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Summary 


An  exact  dynamic  stiffness  element  based  on  higher  order  shear  deformation  theory 
and  extensive  use  of  symbolic  algebra  is  developed  for  the  first  time  to  carry  out  buckling 
analysis  of  composite  plate  assemblies.  The  principle  of  minimum  potential  energy  is 
applied  to  derive  the  governing  differential  equations  and  natural  boundary  conditions. 
Then  by  imposing  the  geometric  boundary  conditions  in  algebraic  form  the  dynamic 
stiffness  matrix,  which  includes  contributions  from  both  stiffness  and  initial  pre-stress 
terms,  is  developed.  The  Wittrick-Williams  algorithm  is  used  as  solution  technique  to 
compute  the  critical  buckling  loads  and  mode  shapes  for  a  range  of  laminated  composite 
plates.  The  effects  of  significant  parameters  such  as  thickness-to-length  ratio,  orthotropy 
ratio,  number  of  layers,  lay-up  and  stacking  sequence  and  boundary  conditions  on  the 
critical  buckling  loads  and  mode  shapes  are  investigated.  The  accuracy  of  the  method 
is  demonstrated  by  comparing  results  whenever  possible  with  those  available  in  the 
literature. 
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1  Introduction 

Several  methodologies  have  been  developed  over  the  years  to  solve  the  elastic  stability 
problem.  A  simplified  approach  to  calculate  the  ith  critical  load,  is  to  consider  the  critical 
load  as  the  load  at  which  more  than  one  infinitesimally  adjacent  equilibrium  configu¬ 
rations  exist  that  can  be  identified  with  the  ith  bifurcation  point  (Euler’s  method)  [1]. 
In  a  linearized  structural  stability  analysis,  the  determination  of  the  critical  load  leads 
to  a  linear  eigenvalues  problem.  The  bifurcation  method  can  be  successfully  used  par¬ 
ticularly  for  plates,  when  the  critical  equilibrium  configuration  shows  a  slight  geometry 
change  as  the  critical  buckling  load  is  reached.  However,  as  explained  by  Leissa  [2], 
linearized  stability  analysis  is  meaningful,  if  and  only  if,  the  initial  in-plane  loading  does 
not  produce  an  out-of-plane  deformation.  Furthermore,  there  are  many  cases  in  which 
Euler’s  method  may  fail,  in  particular  when  thin-walled  structures  like  shells  exhibit 
the  snap-buckling  phenomenon.  In  such  cases,  the  most  general  approach,  based  on  the 
solution  of  the  complete  equilibrium  and  stability  equations  [3,4]  is  preferred. 

Amongst  a  wide  class  of  methodologies  employed  to  analyze  the  elastic  stability  of  ad¬ 
vanced  composite  structures,  the  DSM  is  probably  the  most  accurate  and  computation¬ 
ally  efficient  option.  The  DSM  based  on  Levy-type  closed  form  solution  for  plates  [5]  is 
indeed  an  exact  approach  to  the  solution  procedure.  Wittick  [6]  laid  the  groundwork  of 
the  DSM  for  plates.  The  basic  assumption  in  this  work  is  that  the  deformation  of  any 
component  plate  varies  sinusoidally  in  the  longitudinal  direction.  Using  this  assump¬ 
tion,  a  stiffness  matrix  may  be  derived  that  relates  the  amplitudes  of  the  edge  forces  and 
moments  to  the  corresponding  edge  displacements  and  rotations  for  a  single  component 
plate.  For  the  exact  DSM,  this  stiffness  matrix  is  derived  directly  from  the  equations  of 
equilibrium  that  describe  the  buckling  behavior  of  the  plate.  Essentially,  Wittrick  [6]  de¬ 
veloped  an  exact  stiffness  matrix  for  a  single  isotropic,  long  flat  plate  subject  to  uniform 
axial  compression.  His  analysis  basically  used  classical  plate  theory  (CPT).  Wittrick 
and  Curzon  [7]  later  extended  this  analysis  to  account  for  the  spatial  phase  difference 
between  the  perturbation  forces  and  displacements  which  occur  at  the  edges  of  the  plate 
during  buckling  due  to  the  presence  of  in-plane  shear  loading.  This  phase  difference  was 
accounted  for  by  defining  the  magnitude  of  these  quantities  using  complex  quantities. 
Wittrick  [8]  then  extended  his  analysis  further  to  consider  flat  isotropic  plates  under 
any  general  state  of  stress  that  remains  uniform  in  the  longitudinal  direction  (i.e.,  com¬ 
binations  of  bi-axial  direct  stress  and  in-plane  shear).  A  method  very  similar  to  that 
described  in  [6]  was  also  presented  by  Smith  in  [9]  for  the  bending,  buckling,  and  vibra¬ 
tion  of  plate-beam  structures.  Following  these  developments,  Williams  [10]  presented 
two  computer  programs,  GASVIP  and  VIPAL  to  compute  the  initial  buckling  stress  of 
prismatic  plate  assemblies  subjected  to  uniform  longitudinal  stress  or  uniform  longitu¬ 
dinal  compression,  respectively.  GASVIP  was  used  to  set  up  the  overall  stiffness  matrix 
for  the  structure,  and  VIPAL  demonstrated  the  use  of  substructuring.  Next,  Wittrick 
and  Williams  [11]  reported  on  the  VIPASA  computer  code  for  the  buckling  analysis  of 
prismatic  plate  assemblies.  This  code  allowed  for  analysis  of  isotropic  or  anisotropic 
plates  using  a  general  state  of  stress  (including  in-plane  shear).  The  complex  stiffnesses 
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described  in  [12]  were  incorporated  in  VIPASA,  as  well  as  allowances  were  made  for 
eccentric  connections  between  component  plates.  This  code  also  implemented  an  algo¬ 
rithm,  referred  to  as  the  Wit  trick- Williams  algorithm  [13]  for  determining  any  buckling 
load  for  any  given  wavelength.  The  development  of  this  algorithm  was  necessary  be¬ 
cause  the  complex  stiffnesses  described  above  are  transcendental  functions  of  the  load 
factor  and  half  wavelength  of  the  buckling  modes  of  the  structure  which  make  a  deter¬ 
minant  plot  cumbersome  and  unfeasible.  Viswanathan  and  Tamekuni  [14, 15]  presented 
an  exact  DSM  based  upon  CPT  for  the  elastic  stability  analysis  of  composite  stiffened 
structures  subjected  to  biaxial  in-plane  loads.  The  structure  was  idealized  as  an  as¬ 
semblage  of  laminated  plate  elements  (flat  or  curved)  and  beam  elements.  Tamekuni, 
and  Baker  extended  this  analysis  in  [16]  considering  long  curved  plates  subject  to  any 
general  state  of  stress,  together  with  in-plane  shear  loads.  Anisotropic  material  prop¬ 
erties  were  also  allowed.  This  analysis  utilized  complex  stiffnesses  as  described  in  [12]. 
The  works  described  in  [9, 13, 16]  are  more  or  less  similar.  The  differences  are  discussed 
in  [11].  Williams  and  Anderson  [17]  presented  modifications  to  the  eigenvalue  algorithm 
described  in  [13].  Further  modifications  presented  in  [17]  allowed  the  buckling  mode 
corresponding  to  a  general  loading  to  be  represented  as  a  series  of  sinusoidal  modes  in 
combination  with  Lagrangian  multipliers  to  apply  point  constraints  at  any  location  on 
edges,  the  DS  matrix  for  laminated  composite  plates  for  buckling  analysis.  This  useful 
extension  is  of  considerable  theoretical  and  computational  complexity  as  will  be  shown 
later.  The  research  is  particularly  relevant  when  analysing  thick  composite  plates  for 
their  buckling  characteristics. 
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2  Theoretical  formulation 


2.1  Displacement  field  and  derivation  of  governing  differential  equa¬ 
tions 


In  the  derivation  that  follows,  the  hypotheses  of  straightness  and  normality  of  a  trans¬ 
verse  normal  after  deformation  are  assumed  to  be  no  longer  valid  for  the  displacement 
field  which  is  now  considered  to  be  a  cubic  function  in  the  thickness  coordinate,  and 
hence  the  use  of  higher  order  shear  deformation  theory  (HSDT).  This  development  is  in 
sharp  contrast  to  earlier  developments  based  on  CPT  and  FSDT  and  no  doubt  a  signif¬ 
icant  step  forward.  The  deformation  pattern  through  thickness  of  the  plate  is  shown  in 
Fig.  1.  A  laminated  composite  plate  composed  of  Ni  layers  is  considered  in  order  to  make 
the  theory  sufficiently  general.  The  integer  k  is  used  as  a  superscript  denoting  the  layer 
number  which  starts  from  the  bottom  of  the  plate.  The  kinematics  of  deformation  of  a 
transverse  normal  using  both  first  order  and  higher  order  shear  deformation  are  shown 
in  Fig.  1.  After  imposing  the  transverse  shear  stress  homogeneous  conditions  [18,19]  at 
the  top/bottom  surface  of  the  plate,  the  displacements  field  are  given  below  in  the  usual 
form: 


Figure  1:  Kinematic  descriptions  of  FSDT  and  HSDT  for  a  multilayered  plates. 


u  (x,  y,  z,  t)  =  u0  (x,  y,t)  +  z  4>x  (x,  y,  t)  +  a  z3 

V  (x,  y,  z,  t)  =  Vo  (x,  y,t)  +  z  <\>y  (x,  y,  t)  +  a  z3 
w  (x,  y,  z,  t)  =  W0(x,y,t ) 


(  i  ,  .  dwQ(x,y,t)\ 

U>x  {x,y,t)  + - — -  I 

(  ,  ,  ,x  .  dw0(x,y,t)\ 
{x,y,t)  + - — -  I 


(1) 


where  u,  v,  w  are  the  plate  displacement  components  of  the  displacement  vector, 

T]  =  {  U  V  W  }T  (2) 

ci  =  —  3^  whereas  uo,  vo,  wo  are  the  displacement  components  defined  on  the  plate 
middle  surface  D  in  the  directions  x,  y  and  z.  The  principle  of  minimum  potential  energy 
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is  now  applied.  The  variational  statement  at  multilayer  level  is: 

y]<mfe  =  o  (3) 

k= 1 

where  Tlk  is  the  total  potential  energy  for  the  kth  layer  of  the  composite  plate.  The  first 
variation  can  be  expressed  as: 

8Uk  =  SUk  A  SVk  (4) 

where  5Uk  is  the  virtual  potential  strain  energy,  SVk  is  the  virtual  potential  energy  due 
to  external  loadings,  and  assume  the  following  form: 

5Uk  =  f  f  ( SskT  crk\  dflk  dz ,  SVk  =  f  f  (5efx  dXQ  A  8efy  dyo)  dQk  dz  (5) 

J  J  zk  ^  '  j  J  zk 

the  stresses,  a  and  the  strains,  s  vectors  are  expressed  as  follows: 

&  —  {  ® xx  &yy  Txy  Txz  Tyz  |  |  £xx  £yy  yxy  Axz  A yz  }  (6) 

crXo  and  dyo  denote  the  in-plane  initial  stresses.  The  non-linear  strains  £xlx  and  £yly  are 
approximated  with  the  Von  Karman’s  non-linearity: 

si  =  \(w,xf  4„=\(w,yf  (7) 

The  subscript  T  signifies  an  array  transposition  and  8  the  variational  operator.  Consti¬ 
tutive  and  geometrical  relationships  are  defined  respectively  as: 

ak  =  Ckek  £=T>r)  (8) 

~  k 

where  C  is  the  plane  stress  constitutive  matrix  and  is  the  differential  matrix  (see 
Appendix  A  for  details).  Substituting  Eq.  (42)  into  the  Eq.  (48)  and  imposing  the 
condition  in  Eq.  (46),  the  equations  of  motion  are  obtained  after  extensive  algebraic 
manipulation  as: 

8uo  :  All  Uo,xx  +  Al2  Vo,yx  A  Aiq  (uo  ,yx  T  Vo,xx  )  +  Bn  0  x,xx  A  B 12  4>  y,yx  A  Biq  (( f>x,yx  4"  (fy,xx  )  +  En  C2  <f>  x,xx 
A  E ii  C2  UJo,xxx  +  E12  C2  (j)y,yx  A  E12  C2  Wo,yyx  A  Eiq  C2  4>x,yx  T  Eiq  C2  4>y,xx  A  2  Eiq  C2  VJo,xyx  H-  Aiq  Uo,xy 
+  A26  Vo ,yy  A  Aqq  (uo,yy  A  Vo ,Xy)  A  -016  f>x,xy  +  B2Q  f^y^y  +  Bqq  ( <fx,yy  A  4>y,xy)  A  E 12  C2  {<fx,xy  A  VJo,xxy) 
+  E26  C2  (< fy,yy  A  Wo,yyy )  +  Eqq  C2  ( (fx  ,yy  A  (fy^xy  A  2  KJo ,xyy^)  —  0 

(5^0  •  Aiq  Uq^xx  A  A20  Vo ,yx  A  Aqq  {u>o,yx  A  Vq,xx)  A  Biq  (f)x,xx  A  B20  4>y,yx  A  Bqq  {(f)x,yx  A  (fy^x)  H-  Eiq  C2  (fx,xx 
A  Eiq  C2  VJ 0,xxx  A  E2Q  C2  (j)y,yx  A  -E'26  C 2  VJo,yyx  H“~  Eqq  C2  (fx,yx  A  Eqq  C2  (fy,xx  A  2  Eqq  C2  VJo,xyx  A  A 12  Uo :xy 
A  A22  Vo,yy  A  A2Q  {uo,yy  ~h  Vo,xy )  A  B 12  (fx,xy  4"  B 22  8y ,yy  4"  B2Q  ( 4>x,yy  A  (fry^xy')  A  E 12  C2  ((f>x,xy  A  U)o ,xxy^) 
A  E22  C2  ((fy,yy  A  VJQ^yyy)  A  -E'26  C2  (< fix,yy  A  (fy  ,xy  4"  2  VJo,xyy )  —  0 


(9) 
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Swo  :  A44  fay,y  E  wo,yy)  E  A45  fax,y  E  VJq ,xy)  E  B44  Cl  ( 4>y,y  E  Wo,yy)  +  -D45  Cl  (( f>x,y  E  VJo ,xy) 

E  A45  fay,x  E  ^0,X2/)  +  ^55  fax,x  E  VJq,xx )  +  D45  Cl  (</>y,cc  +  ^0,ccy)  +  B55  Cl  ( (fx,x  E  VJq,XX) 

+  D 44  Ci  ((fy,y  E  ,yy)  E  B45  Cl  ((j)Xjy  H“~  ,xy)  E  F44  Ci  ((fy,y  H-  ,yy)  E  B45  C\  ((J)Xjy  E  VJo ,ccy) 

+  B45  Cl  fayjX  H-  ^0 ,ccy)  +  B55  Cl  (( f>x,x  E  VJo,xx^)  H-  B45  Ci  fay,x  H-  ^0 ,ccy)  H-  B55  Ci  fax,x  E  VJo,xx^) 

Bn  C2  Uo,xxx  -E'12  C2  VojXXy  -E/16  C2  fao,xxy  H-  'C0,cccccc)  Bn  C2  4*x,xxx  B12  C2  8y,xxy 
Bi6  C2  (< fx,xxy  E  (fry, xxx  )  -  Hu  cl  (4>x  ,xxx  E  VJq,xxxx)  H\2  C2  fax  ,xxy  E  V)o,xxyy') 

-f^l6  C2  (< fix,xxy  E  (fy,xxx  E  2lCo ,xxxy)  2  Bi6  C2  Uo,XXy  2  B26  C2  Vo,xyy  2  Eqq  C2  fao,xyy  E  'C0,£cccj/) 

2  -Fl6  C2  <fx,xxy  2  -F26  C2  4>y,xyy  2  B36  C2  fax,xyy  H-  (fy,xxy^)  2  i^l6  C2  fax,xxy  E  VJo,xxxy^) 

2  H26  C2  fay,xyy  E  ,xyyy)  2  Hqq  C2  ((fx ,xyy  E  (j)y,xxy  E  2  VJo  ,xxyy)  E42  C2  ^0 ,xyy  E22  C2  ^0 ,yyy 
E2Q  C2  ( U0,yyy  E  Vo^xyy')  B12  C2  (fx,xyy  F22  C2  f^y^yyy  B26  C2  (< fix,yyy  E  (fy,xyy^) 

H\2  C2  fax,xyy  E  VJo ,xxyy^)  H22  C2  fay,yyy  E  VJo ,yyyy^)  2  B26  C2  fax,yyy  E  8y-,xyy  H-  2  VJo^xyyy') 

—  JV*0  VJq,xx  E  -/Vj/o  ^0,yy 

S(f)x  •  Bn  uojXX  E  B\2  vq ,yX  E  Biq  (u o  ,£/cc  +  'Co,cccc  )  +  £>n  <t>  x,xx  E  D 12  4>y,xy  E  E)\q  (( fix,yx  E  4>y,xx) 

E  F\l  C2  ( 4>x,xx  E  ^0,cccccc)  +  B12  C2  ( <fty,yx  E  Wo ,yyX)  E  F±Q  C2  (0x,ycc  H-  4*y,xx  H-  2lCo,ccycc) 

“b  Bi6  UojXy  E  B 26  V0,yy  E  Bqq  [liQ^yy  E  Vq ,xy)  E  D±q  (f)x,Xy  E  E) 26  8y,yy  E  E)qq  ( (fx,yy  E  (fy,xy^) 

E  F\q  C2  ( 4>x,xy  E  VJo ,XXy^)  H-  F2Q  C2  {^(fy^yy  E  'tC0,j/j/y)  H"  -^66  C2  (( f>x,yy  E  8y,xy  E  2  ICo ,Xyy^) 

E  En  C2  UojXX  H-  B12  C2  'Co,yx  +  Bi6  C2  ( U0,yX  E  'Co,xcc)  H-  Bn  C2  (f>x,xx  H-  B12  C2  4>y,xy  E  Fiq  C2  {(fx,yx  E  8y,xx^) 
E  Hu  cl  (fa  ,xx  E  Wo:XXX)  E  H12  C2  (4>y,yx  E  Wo:yyx)  +  i^l6  C2  (0x  H-  (fy^x  E  2ico,xyx) 

+  Bi6  C2  ^0,cc?/  +  E26  C2  +  -^66  C2  (^0,yj/  +  ^0,a;y)  +  Bl6  C2  (fx,xy  +  B26  C2  <fy,yy  E  Fqq  C2  (4>x,yy  E  4>y ,xy) 

E  EIiq  C2  (( f>x,xy  E  Wq ,xxy^)  H-  B26  C2  (8y,yy  E  ^0 ,yyy^)  E  Hqq  C2  (c f>x,yy  E  8y,xy  E  2  WojXyy^j 
4I45  (c f>y  E  2  4Co,y)  4I55  (cfx  E  2  ico,x)  2  E) 45  ci  ~b  2  ico,y)  2  E) 55  ci  H-  2  ico,cc) 

-  B45  Cl  (fa  E  2  tco,y)  -  B55  Cl  (fa  E  2  ic0,x)  =  0 

•  Bi6  U0,XX  E  B2Q  Vo ,yX  E  Bqq  (lLo,yx  E  Vq,xx )  H-  Bl6  (f>x,xx  E  D2Q  (fy,xy  E  E)qq  fax,yx  E  (fy,xx ) 

+  Bi6  C2  (( fix, XX  E  VJojXXX^  E  F2Q  C2  (4>y,yx  E  VJQ,yyx)  E  Fqq  C2  ((fx,yx  E  (fy^xx  E  2  WQ,xyx) 

E  B12  V>o,xy  E  B22  Vq ,yy  E  B2Q  (v>o,yy  E  Vo ,xy^)  E  E)  12  (fx,xy  E  E) 22  f^y,yy  H-  B2Q  (8x,yy  E  (fy,xy') 

E  F12  C2  (4>x,xy  E  VJo ,XXy^)  H-  E22  C2  fay^y  E  VJo,yyrf)  E  F2Q  C2  ((fx,yy  E  (fy-tXy  E  2  ICf) ,xyy^) 

E  Eiq  C2  UojXx  H-  E2Q  C2  Vo,yx  E  Eqq  C2  (uo,yx  H-  'C0,xcc)  E  Fiq  C2  4>x,xx  H-  F2Q  C2  (fy,xy  E  Fqq  C2  (<fx,yx  E  4>y,xx ) 
E  Hiq  C2  (4>x,xx  E  VJo^xxx )  E  H2Q  C2  (<fy,yx  H-  V)Q,yyx)  E  Hqq  C2  fax  ,yx  H-  (fy,xx  E  2lCo,xyx) 

+  B12  C2  Uo,xy  E  E22  C2  Vq ,yy  E  E2Q  C2  fao,yy  E  Vq ,Xy)  E  F12  C2  4>x,xy  +  F22  C2  4>y,yy  E  F2Q  C2  fax,yy  E  fa ,xy) 
E  H12  C2  fax,xy  E  VJo ,xxy)  E  H22  C2  fay,yy  E  VJo,yyy^}  H-  H2Q  C2  fax,yy  E  4>y,xy  E  2  Wo,xyy ) 

-A44  ((fy  E  2  4Co,y)  4^45  ((fx  H-  2  lCo,cc)  2  F) 44  Cl  ((fy  E  2  lCo,y)  2  F) 45  Ci  +  2  lCo,x) 

—  F44  Ci  (0y  +  2  ico,y)  —  B45  Ci  (0X  +  2  icq,®)  =  0 
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The  natural  boundary  conditions  are: 

SlLo  •  Nxx  —  All  Vq,x  +  Bn  (f>x,x  A  En  C2  (fx,x  A  Ell  C2  Wo,xx  H-  -^4.12  Vo ,y  A  £l2  (fy,y  A  E12  C2  H-  E12  C2  WQ,yy 
A  Aiq  Uo,y  A  Aiq  Vq,x  A  Biq  <$>x,y  A  B\q  (j)y,x  A  E\q  C2  (fx,y  A  Eiq  C2  +  2  Eiq  C2  Wo,xy 

Sv 0  :  Mxy  —  Aiq  Uo,x  A  Biq  (f)x,x  A  Eiq  C2  (frx,x  A  Eiq  C2  Wo, XX  +  A2Q  Vo ,y  A  B2Q  <j)y,y  A  E20  C2  4>y,y  A  E20  C2  Wo,yy 

A  Aqq  Uo,y  A  Aqq  Vo,x  A  Bqq  4>x,y  A  Eqq  C2  4>y,x  +  Eqq  C2  4>x,y  A  Eqq  C2  (j)y,x  A  2  Eqq  C2  Wo,xy 


Swo  : 


Qx  m  Hu  C2  4>  x,xx  A  Hu  C2  Wo, xxx  T  Ell  C2  Uo,xx  A  Fn  C2  4>x,xx  T  E12  C2  Vo,yx  A  E12  C2  4>y,yx 

A  Hi 2  C2  4>y,yx  +  Hi 2  C2  Wo,yyx  A  2  £16  C2  Uo,xy  A  2  Fiq  C2  <t>x,xy  A  2  H1QC2  4>x,xy  +  £l6  C2  UQ,yx 
A  Eiq  C2  vo, xx  A  Fiq  C2  4>  x,yx  +  h16  cl  4>  x,yx  +  Hie  cl  4>y  ,xx  H-  2  Hiq  C2  Wo,xxy  A  2  £/26  C2  VQ,yy 
A  2  £26  C2  (fy^yy  A  2  H2QC2  Wo,yyy  A  4  HqqC2  Wo,xyy  H-  2  H2Q  C2  (fx,yy  A  2  H2Q  C2  (j)y,xy  A  2  Eqq  C2  Uo ,yy 
A  2  Eqq  C2  Vo,xy  ~"b  2  Fqq  C2  (fx,yy  A  2  Fqq  C2  (fy^xy  2  £>45  Ci  (j)y  2  £>45  Cl  Wo,y  £45  Cl  (fy 
—  £45  Ci  —  -4.55  4>x  ~  ^-55  ~  Eqq  Cl  (f>x  ~  2  Cl  Ido, x  “  £55  C?  (j)x  ~  Fqq  c\  Wq,x 


S(f)y  •  MiXx  —  £>n  (fx,x  +  ££l.l  C2  (fx,x  A  H Ll  C2  Wo, xx  A  Bn  Uo,x  A  En  C2  T  2  £ll  C2  4>x,x  A  Fn  C2  Wo, xx 

A  £11  C2  +  £>12  ^0,2/  +  £>12  4>y,y  A  £12  C2  (j)y,y  A  £12  C2  Wo,j/2/  +  £l2  C2  Vo,y  A  £l2  C2  4>y,y  A  H\2  (\)y,y 

A  H L2  C2  Wo,yy  +  £?16  +  £>16  'CO,^  +  £>16  +  £>16  0j/,cc  +  £l6  c2  +  £l6  c2  4>y,x  A  2  £16  C2 

+  £7ie  C2  'W'O,?/  +  Eiq  C2  fO,x  +  £16  C2  4>x,y  A  Fiq  C2  0y,cc  +  ££l_6  C2  (fx,y  A  H\q  C2  4>y,x  +  2  ££.6  C2  WQ,xy 


S(f)x  •  Adxy  —  £>16  4*x,x  A  Hiq  C2  4>x,x  H-  Hiq  C2  Wo, xx  H-  -£l6  Vo,x  A  Eiq  C2  x  A  2  £16  C2  (fx,x  A  Fiq  C2  Wo, xx 

A  Fiq  C2  Wo, xx  A  £>26  Vo,y  A  £>12  <fy,y  A  £26  C2  4>y,y  A  £26  C2  Wo,yy  A  E20  C2  Vo ,y  A  £26  C2  4>y,y  A  H2Q  C2  4>y,y 
A  H2Q  C2  Wo,yy  A  Bqq  Uo,y  A  Bqq  Vq,x  A  £>66  (fx,y  A  Dqq  (])y,x  A  Fqq  C2  (fx,y  A  Fqq  C2  A  2  £66  C2  Wo,xy 
A  Eqq  C2  V>o,y  A  Eqq  C2  Vq,x  A  Fqq  C2  (fx,y  A  Fqq  C2  (j)y,x  A  Hqq  C2  (fx,y  A  Hqq  C2  (fy,x  A  2  Hqq  C2  Wo,xy 


Swo  ,X  •  £ CCX  -  i£ll  C2  (fx,x  A  Hu  C2  lCo,ccx  +  £11  C2  ltO,cc  +  £11  C2  (fx,x  +  £12  C2  Vo ,y  A  £12  C2  4>y,y  A  ££2  C2  4>y,y 

A  H12  C2  Wo,yy  A  £l6  C2  Uo,y  A  £l6  C2  fOjX  +  £l6  C2  4>x,y  A  Fiq  C2  <fty,x  A  H\q  C2  4>x,y  A  Hiq  C2  (fy,x 
A  2  Hiq  C2  Wo,xy 

(10) 


where  the  suffix  after  the  comma  denotes  the  partial  derivative  with  respect  to  that 
variable  and 


(Aij ,  Bjj ,  Dij ,  £/2j ,  F{j ,  HXJ  j 


(h,  h,h,  h,  h,  h) 


(11) 


are  laminate  stiffnesses  and  rotatory  inertial  terms,  respectively  with  i  and  j  varying 
form  1  to  6.  The  in-plane  loadings  can  be  defined  as  NXo  =  A  NXo  and  Nyo  =  A  Nyo , 
where  NXo,  Nyo  are  the  initial  in- plane  loadings  and  A  is  a  scalar  load  factor,  ci  has 
already  been  defined  (see  Eq.  (37))  and  C2  =  —  A. 
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2.2  Dynamic  stiffness  method 

Once  the  equations  of  motion  and  the  natural  boundary  conditions,  i.e.,  Eqs.  (9)  and 
(10)  above  are  obtained,  the  classical  method  to  carry  out  an  exact  buckling  analysis 
of  a  plate  consists  of  (i)  solving  the  system  of  differential  equations  in  Navier  or  Levy- 
type  closed  form  in  an  exact  manner,  (ii)  applying  particular  boundary  conditions  on 
the  edges  and  finally  (iii)  obtaining  the  stability  equation  by  eliminating  the  integration 
constants  [20-23].  This  method,  although  extremely  useful  for  analysing  an  individual 
plate,  it  lacks  generality  and  cannot  be  easily  applied  to  complex  structures  assembled 
from  plates  for  which  researchers  usually  resort  to  approximate  methods  such  as  the 
FEM.  In  this  respect,  the  dynamic  stiffness  method  (DSM),  which  is,  in  many  ways, 
analogous  to  FEM  has  no  such  limitations  and  importantly  it  always  retains  the  exactness 
of  the  solution  even  when  applied  to  complex  structures.  This  is  because  once  the 
dynamic  stiffness  matrix  of  a  structural  element  is  obtained  from  the  exact  solution  of 
the  governing  differential  equations  and  it  can  be  offset  and/or  rotated  and  assembled  in 
a  global  DS  matrix  in  the  same  way  as  the  FEM.  This  global  DS  matrix  thus  contains 
implicitly  all  the  exact  critical  buckling  loads  of  the  structure  which  can  be  computed 
by  using  the  well  established  algorithm  of  Wit  trick- Williams  [13]. 

A  general  procedure  to  develop  the  dynamic  stiffness  matrix  of  a  structural  element  is 
generally  summarized  as  follows: 

(i)  Seek  a  closed  form  analytical  solution  of  the  governing  differential  equations  of  the 
structural  element. 

(ii)  Apply  a  number  of  general  boundary  conditions  in  algebraic  forms  that  are  equal 
to  twice  the  number  of  integration  constants;  these  are  usually  nodal  displacements 
and  forces. 

(iii)  Eliminate  the  integration  constants  by  relating  the  amplitudes  of  the  harmonically 
varying  nodal  forces  to  those  of  the  corresponding  displacements  which  essentially 
generates  the  dynamic  stiffness  matrix,  providing  the  force-displacement  relation¬ 
ship  at  the  nodes  of  the  structural  element. 

Referring  to  the  equations  of  motion  Eqs.  (9),  an  exact  solution  can  be  found  in  Levy’s 
form  for  symmetric,  cross  ply  laminates.  For  such  laminates  B  —  E  —  0,  and  C\§  = 
=  C45  =  0  and  the  out-of-plane  displacements  are  uncoupled  from  the  in-plane  ones. 

2.3  Levy-type  closed  form  exact  solution  and  DS  formulation 

The  solution  of  Eqs.  (9)  related  to  the  out-of-plane  displacements  is  sought  as: 

00  00 

w°(x,y,t )  =  y  Wm{x)elult  sm(uy),  cj)x(x,y,t)  =  y  $xm(x)  elut  sin(ay), 

ra= 1  ra= 1 

00  '  ' 

<f>y(x,y,t)  =  y  ®ym(x)  elut  cos  (ay) 

m= 1 
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where  uj  is  the  unknown  circular  or  angular  frequency,  a  =  and  m  =  1,2,. . . ,  oo. 
Equation  (12)  is  the  so-called  Levy’s  solution  which  assumes  that  two  opposite  sides  of 
the  plate  are  simply  supported  (S-S),  i.e.  w  =  (j)x  =  0  at  y  =  0  and  y  —  L.  Substituting 
Eq.  (12)  into  Eqs.  (9)  a  set  of  three  ordinary  differential  equations  is  derived  which  can 
be  written  in  matrix  form  as  follows: 


£11  £12  £13 

■  Wm  ' 

"  0  " 

£21  £22  £23 

3>* 

= 

0 

£31  £32  £33 

0 

where  Cij  —  1,  2,  3)  are  differential  operators  and  given  by: 

£n  =  —CiV^Hn  +  V2  (A55  +  2C2H55  4-  C2E55  +  2a2 c\H  12  +  4o2 c2Hqq  +  AA^o)  —  <a2  ^44  +  2C2D44 
+  C2-E44  +  cf2  c2H 22  +  AA^yo^ 

£12  =  Px  (—ci Eh  ~  c2#n)  +  Vx  (A55  +  2C2H55  +  o2ciF12  +  C2E55  +  2a2ciFQQ  +  a2  c2H  12  +  2o2c2i?66) 
+  (—ci  Fn  —  c2  Hu)  T>\ 

£l3  =  —  OL  (A44  +  02(21)44  +  C2E44)  +  02Ci  (#22  4“  C1H22))  +  D2  (o'CiEi2  4"  2qlC\Fqq  +  OC2Hi2  +  2oc2#66) 
£21  =  ci#8(#n  +  ciUn)  +  Vx  (— 4L55  —  <^2(21)55  +  C2E55)  —  o2ci(F12  +  2#66  +  ci#i2  +  2ci#66)) 

£22  =  —^55  —  <^2(21)55  +  C2E55)  +  D2  (Dn  +  2ciDn  +  c2Hn)  —  a2  (D66  +  2ciD66  +  c2!^) 

£23  =  Px  (— 0D12  —  cvDqq  —  20C1D12  —  2ociD66  —  ac2Hi2  —  nc2#66) 

£31  =  —  (A44  +  C2 (2D44  +  C2D44)  +  02Ci(D22  +  C1U22))  +  aCiD2(Di2  +  2#66  +  C\H\2  +  %CiHqq) 

£32  =  (FDx(Di2  +  D66  +  ci(2Di2  +  2#66  +  ci(#i2  +  #66))) 

£33  =  —  AL44  —  C2  (2D44  +  C2D44)  —  02(D22  +  Ci(2D22  +  Cl  #22))  +  #2(#66  +  Ci(2D66  +  Cl  #66)) 

(14) 

where  and  #q-,  C^,  Dq-,  #^-,  #q-,  #q  have  already  been  defined  in  Eq. 

(57).  Expanding  the  determinant  of  the  matrix  in  Eq.  (13)  the  following  differential 
equation  is  obtained: 


(0-1  D8  +  0-2  Px  +  0-3  D4  4-  0-4  +  0,5)  4/  —  0  (15) 

where 

«  =  (16) 

Using  a  trial  solution  eA  in  Eq.  (63)  yields  the  following  auxiliary  equation: 

0-1  A8  H-  0/2  A6  +  0-3  A4  4-  0-4  A^  4-  0-5  =  0  (17) 

Substituting  p  —  A2,  the  order  polynomial  of  Eq.  (17)  can  be  reduced  to  a  quartic 
as: 

oq  p^  4~  oq  4“  0-3  P'2  +  0-4  4~  0-5  =  0  (18) 
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the  four  roots  for  the  quartic  equation  are  given  by: 


Mi  =  sx-\ 
1 

M2  —  —si  +  2 

M3  =  -sx-\ 

1 

£4  —  —  S\  +  - 


i 

—  S  5  +  S2  — 

58 

56 

4  y^9 

3  a\  57 

1 

—  S5  +  S2  — 

58 

56 

4  y/SQ 

3  a\  57 

/  ~~ 
—  S5  +  S2  + 

58 

56 

4  y^9 

3  a\  57 

—  55  +  52  + 


58 


<§6 


4^/59  3ai  57 


+  \y/*» 
+  \^ 


(19) 


where 


a\  4  a2  22  2 

51  =  - — ,  52  =  —  — o  +  — o,  53  =  2a3  -  9a2  a3  a4  -  72 a2  a5  +  27 a2  a5  +  27 a\  a4, 
4a2  3af  2af 


zo  l/^53  +  a/53  —  4  54^  3 

54  =  a3  —  3  <22  <24  +  12  a  1  <25,  55  =  — 

oq 


56 


=  .s 4,  s7  =  \/323 


«5  Ml,  «8  = 


32 

<22  \  3  4  a2  03  8  04 

a\  J  a?  ai  ’ 


_  «2  s@ 

S9  —  S5  +  —  +  - - 

2  o  57  ai 

(20) 


The  explicit  form  of  the  polynomials  coefficients  aj  ( j  =  1,2,  3,4,  5)  are  given  in  Ap¬ 
pendix  B.  Some  pair  or  pairs  of  complex  roots  may  occur  when  computing  /ij  (j  =  1,  2,  3, 4), 
but  the  amplitude  of  the  displacements  Wm  (x) ,  4>Xrn  (x) ,  <I>yrn  (x)  are  all  real,  whilst  the 
associated  coefficients  can  be  complex.  As  complex  roots  always  occur  in  conjugate  pairs, 
the  associated  coefficients  will  also  occur  as  conjugates.  The  solution  of  the  system  of 
ordinary  differential  equations  in  Eq.  (13)  can  thus  be  written  as: 

Wm  (x)  =  A\  e+w  *  +  A2  e_w  x  +  A3  e+^2  *  +  A4  e~^  x 
+  A5  e+M3  x  +  Ag  e_/43  x  +  A7  e+M4  x  +  A8  e_/44  x 


(x)  =  Bi  e+/41  x  +  B2  e “/il  *  +  £3  e+/42  x  +  B4  e“ ^2  x' 
+  B5  e+M3  x  +  B6  e_M3  x  +  B7  e+/44  x  +  £8  e_/44  x 


(21) 


^  (®)  =  Ci  e+M1  x  +  C2  e_M1  *  +  (73  e+M2  x  +  C4  e“w  x 
+  C75  e+M3  X  +  C6  e_M3  *  +  C7  e+/i4  *  +  C8  e"^4  * 

where  Ai  —  A8,  B\  —  B%,  C\  —  C8,  are  three  sets  of  integration  constants.  The  sets  of 
constants  are  not  all  independent.  Only  one  set  of  eight  constants  are  needed  to  relate 
each  set.  Constants  B\  —  B$  are  chosen  to  be  the  independent  base.  By  substituting 
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Eqs.  (66)  into  (13)  the  following  relationships  are  obtained  using  symbolic  computation: 


A i  =  Si  B i, 

A2  =  —8 1  B), 

Ci 

=  7i  £i> 

C2  = 

-71  £ 

As  =  82  Bs, 

A4  =  —82  b4, 

c3 

=  72  £3, 

c4  = 

-72  £ 

As  =  £3  Bs, 

Aq  =  —8s  Be, 

C5 

=  73  £5, 

C6  = 

-73  £< 

A7  =  84  B7 , 

Ag  =  —84  B8 , 

C7 

cq 

II 

C8  = 

-74  Bt 

where 

Si  =  —  [  —  A55a2  £>22  —  2A55C2-D44 

—  2a2 c2  £>22  £55  _  4c2£>44£>55 

—  a4  £>22  £>66  —  2a2 

c2  £44  Dqq  — 

2  A55a2( 

—  4c*2  ci  C2  £>55  ^22  —  2a4ci£>66F22  —  A55c|f44  —  2c2D55F44  —  oc2c2DqqF 44  —  a2  c2D22F55  —  2c2£44F55 

—  2a2 ci c2 F22 F55  ~  c2 £ 44 F 55  —  2a4ci£>22£166  —  4a2ciC2£44£66  —  4a4c2F22£66  —  2a2cic2F44F66 


(22) 


—  2o.Ac\f22Hqq  —  a2c2C2 £44-^66  —  cd ci H 22 H 55  —  A44  (A55  +  2c2D55  +  c2F55  +  a2(£>66  +  2ciF66  +  c2Hqq)^ 

+  A44 (£>n  +  ci(2Fn  +  ciifu))^2  +  (A55  (£>66  +  ci(2F66  +  c\Hqq))  +  2c2(£>n£>44  +  D55Dqq  +  ci(2£>44Fh 

+  2D55F66  +  Ci  £>44 +  Ci  £>55/^66))  +  C2(£>ll£44  +  DqqFqs  +  Ci(2£hF44  +  2F55F66  +  Ci  F44F1  1  +  C1F55F66)) 

—  a2(£22  -  £>h(£>22  +  cl  (2F22  +  C1F22))  +  2£i2(£66  +  ci(2Fi2  +  2F66  +  ci(Fi2  +  H66)))  +  ci(4Fi2(£>66  +  ciFi2) 

—  £>22(2Fh  +  ci  i£n)  +  ci(8Fi2F66  +  2DqqH12  —  2Fn(2F22  +  c\H22)  +  ci(  —  2F22Hn  +  Fi2(4(£i2  +  F66)  +  ^.^12) 


—  C1F11F22  +  4Fi2i£66  +  2c1Hi2H66))))^ /j,2  —  (Du  +  ci(2Fn  +  c1H11))(Dqq  +  ci(2F66  +  cii£66))^4] 

/  [«2(Fi2  +  £>66  +  ci(2Fia  +  2F66  +  ci  (H\2  +  if66)))/Xi  (a44  +  2c2F44  +  c2F44  +  a2ci(F22  +  C1F22) 

—  ci  (F12  +  2F66  +  ci  if  12  +  2ci  Hqq)/j,2 )  —  m  (A55  +  2c2D55  +  c2F55  +  a2ci  (Fi2  +  2F66  +  c\H12  +  2ci£f66) 

—  ci  (F 1  1  +  ciifu)^2)  (A44  +  c2  (2F44  +  C2F44)  +  ex2  (D22  +  2ci  F22  +  c2H22  )  —  (F66  +  2ciF66  +  c2F66)Mi)] 


7i  =  [g  (a44  A55  +  2A55C2F44  +  2A44C2F55  +  4C2F44F55  +  A44a2£66  +  2a2 c2 £>44 £>66  +  A55a2ciF22  +  2a2cic2£s5£22 
+  a4ci  DqqF22  +  A55C2F44  +  2C2F55F44  +  a2c2F66F44  +  A44C2F55  +  2c2£44F55  +  a2cic2F22£55  +  ^£44  £55 
+  2A44a2ci  Fge  +  4a2ciC2£>44-£66  +  2a4c2F22£66  +  2a2ci  c^F^Fqq  +  A^^cx2  c2H22  +  2a.2c2c2Dr)^H22  +  cx.Ac\DqqH22 
+  a2c2C2F55if22  +  2cxAc\FqqH22  +  A44a2c2if66  +  2  a2c2c2  F44if66  +  a4cf  F22H66  +  a2  c^c^F^Hqq  +  a4c4  H22Hqq 
+  (^55(^12  +  F66)  -  A44(£>ii  +  ci(2Fn  +  ciFu))  -  2c2  (£>11  £>44  -  £)55(£)12  +  £>66)  +  ci(2£>44Fh  -  D55F12  +  ci£44Fh 
+  ci  D55  Hqq))  +  c2(  — F11F44  +  (£12  +  Dqq)F55  —  ci  (2Fn  F44  —  F12F55  +  C1F44F11  +  ciF55H66))  +  ci(A55(Fi2  —  C1F66) 

+  a.2  (  — Fn  F22  +  £>12  (£12  +  2F66  +  ci  (H\2  +  2Hqq))  +  ci(2F22  -  DnH22  -  2F11(F22  +  ciH22)  +  £l2(4F66  +  3ciFi2 
+  4ci  Hqq)  +  ci  (F12  (2F66  +  C1F12)  -  £11  (£22  +  C1F22)  +  2ci  Hi2Hqq))))^  fx2  +  c1(D11(F12  +  2F66)  -  £>12  (£11  +  ci  Jfn) 

-  £>66 (£ll  +  C4  F 1  1  )  +  ci  (£>11  (F12  +  2 F66)  +  Ci  Fn  (-Fi2  +  ciF66)  +  Fn(2F66  +  ci  (Fi2  +  3£f66))))^4)] 

/  \jx,i  (A44A55  —  A44a2£>i2  +  A55CX2 D22  +  2 A55 c2 £>44  —  2a2C2£>i2£>44  +  2 A44C2 £>55  +  2a2C2£>22£>55  +  4c^ £>44 £>55 

22  24  2  24  2 

—  A44CK  £>66  —  2a  c2 £>44 Dqq  —  A44a  C1F12  +  oc  c\D22F\2  —  2a  ciC2£>44£l2  +  2Ag5a  C1F22  —  ot  c\Di2F22  +  4a  cic2Dr,r,F22 


4  222  3  22  22  222 

—  a  c\DqqF22 A^c2F^  —  oc  c2£>12£44  +  2c2£>55£44  —  «  c2£>66£44  —  cn  ci c2Fi2£44  +  A^c2Fr)z)  +  a  c2D22F^r) 

+  2c2£>44F55  +  2a2cic2F22£55  +  c2F44F55  +  2at4ciD22F6Q  +  2a4  c2F22Fqq  +  a4c2 F22F12  +  a4cfF22Fi2  +  A55a2c2F2  2 

—  a4c2 £>i2F22  +  2a2c2c2  £>55^22  —  «4c2F66F22  —  a4cf  Fi2F22  +  oc2  c2c2F55H22  +  A44a2c2F66  +  2a4c2D22H66 
+  2a2c2c2£>44£66  +  3a4  c^F2  2^66  +  Q:2c2c2F44F66  +  a4c4F22F66  —  (a5s(£>66  +  ci(2F66  +  ciF66))  +  c2(D66F5  5 
+  ci(FhF44  +  2F55F66  +  C1F44F11  +  ciF55F66))  +  2c2(ci£>44(Fh  +  ci^n)  +  D55(D66  +  ci(2F66  +  ciF66))) 

+  ci(A44(Fii  +  ciFu)  +  a2  (  —  £12  (£>12  +  2ciFi2)  —  2£>i2£66  +  £)22(£ll  +  ci^n)  +  ci(— 4Fi2£66  +  £ll(2£22  +  ciH22) 

—  D12(H12  +  2Hqq)  +  c1(2F22H11  —  H12(3F12  +  2F66  +  ci^i2)  +  ci  Fi  1  F22  —  4Fi2F66  —  2ciFi2F66)))))  id 
+  ci(Fn  +  ci  1 )  (Dqq  +  ci  (2Fqq  +  ciF66))Mi)] 

(23) 

with  i  =  1,2,  3, 4.  The  procedure  leading  to  Eqs.  (22)  and  (23)  must  be  undertaken 
with  sufficient  care,  because  if  wrong  equations  are  chosen  from  Eq.  (66)  to  obtain  the 
relationship  connecting  different  sets  of  constant,  numerical  instability  can  occur.  When 
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Eqs.  (22)  are  substituted  into  Eqs.  (66)  a  solution  in  terms  of  only  8  constants  can  be 
formulated  for  Wm  (x),  &xm  ( x )  and  4>Vm  (x),  respectively.  Thus 

Wm  (x)  =  Bi  e+ ^  x  -B2Si  e“ ^  x  +  B382  e+^2  x  -  B452  e“M2  x 
+  B5  53  e+w  x  -B683  e“^3  X  +  B784  e+/44  x  -  B8S4  e“M4  x 

4>Xm  (x)  =  b4  e+w  x  +  B2  e_/41  x  +  B3  e+^2  x  +  B4  e“^2  x 

+  B5  e+^3X  +  B6  e~^x  +  B7  e+^x  +  B8  e~^x  ^  ' 


<hym  (x)  =  B4  7l  e+Ml  x-B2  7l  e_Ml  x  +  B3l2  e+^2  x  -  B4l2  e~^  x 
+  B5  73  e+M3  x-B6  73  e_M3  x  +  B7  74  e+^4 *  -  S8  74  e“w  31 

The  expressions  for  forces  and  moments  can  also  be  found  in  the  same  way  by  substi¬ 
tuting  Eqs.  (69)  into  Eqs.  (10)  and  using  symbolic  computation.  In  this  way 

Qx  (x,  y )  —  x  (By  +  B2  e  2  x)(A 55  +  T55  Si  mi  +  2  C2  (E55  +  E55  Mi)  +  c2  (E55  +  E55  Mi) 

+  ci  (a 71  (F12  +  2  E36  +  ci  H12  +  2  ci  1^66)  Mi  —  Mi  (Eli  +  ci  Hu  +  ci  Si  Hu  Mi)  +  (2  E36 

+  2  ci  Hqq  +  ci  Si  H12  Mi  +  4  ci  Si  Hqq  mi)))  + 

eP2  x(b3  +  S4  e  2^2  "XA55  +  A55  ^2  M2  +  2  C2  (E55  +  E55  £2  M2)  +  c2  (E55  +  S2  Fqq  M2) 

+  ci  (a  72  (Ei 2  +  2  Fqq  +  ci  H12  +  2  ci  Hqq)  M2  —  M2  (Eli  +  ci  Hu  +  ci  £2  -Hn  M2)  +  «2  (2  E36 
+  2  Cl  Hqq  +  Cl  S2  H\2  M2  +  4  Ci  <h  Hqq  /l 2)))  + 

eP3  x  4-  56e  2  X)(A55  +  A55  (^3  M3  +  2  C2  (E55  +  E55  S3  M3)  +  C2  (E55  +  S3  Fqq  /I3) 

+  Cl  ( a  73  (El2  +  2  ^66  +  Cl  H\2  +  2  Cl  Hqq)  M3  ~  M 3  (Eli  +  Ci  Sn  +  Cl  S3  Hu  M3)  +  CK2  (2  Fqq 
+  2  Cl  Hqq  +  Cl  S3  H\2  M3  +  4  Ci  S3  Hqq  M3)))  + 

eM4  (E3  +  S4  e  )  (-A55  +  Aqq  S4  M4  +  2  C2  (E55  +  E55  S4  M4)  +  C2  (E55  +  S4  Fqq  M4) 

+  Cl  (a  74  (El2  +  2  ^66  +  Cl  H\2  +  2  Cl  i^66)  M4  —  M4  (Ell  +  Cl  H\\  +  Cl  S4  Hu  M 4)  +  «2  (2  Fqq 
+  2  Cl  #66  +  Cl  54  #12  ^4  +  4 Cl  54  #66  Mi))))  sin  (a  y)  =  Qx  sin  (ay) 


M**  (x,  y)  =  (e'41  *(#1  +  #2  e-2M1  x)(a2  ci  <5i  (#12  +  ci  #12)  +  a 71  (£>12  +  a  (2  #2  +  a  #12))  -  n  1  (#11 
+  ci  (2  Eh  +  ci  Hu  +  <h  Eli  Mi  +  ci  ^1  Hu  Mi)))  + 

e^2  "(S3  +  S4  e-2  ^2  ")(a2  ci  (Ei2  +  ci  Hy2)  +  a  72  (Si2  +  ci  (2  E12  +  ci  Si2))  -  M2  (^11 
+  ci  (2  Eli  +  ci  Ei  +  S2  Eli  M2  +  ci  ^2  Ei  M2)))  + 

e^3  +  S5  e  2  ")(a ?  ci  (^3  (El2  +  ci  S12)  +  ol  73  (-D12  +  ci  (2  El 2  +  ci  S12))  —  M 3  (Ei 

+  ci  (2  Eh  +  ci  Si  1  +  £3  Eh  M 3  +  ci  $3  Sn  ms)))# 

eM4  "  (S7  +  S8  e“2  ^4  ")(a2  ci  S4  (E12  +  ci  H12)  +  a  74  (£>12  +  ci  (2  Ei2  +  ci  H12))  -  M4  (Eh 
+  ci  (2  En  +  ci  Sn  +  S4  En  M4  +  ci  £4  S11  M4))))  sin  (a  y)  =  Mxx  sin  (a  y) 
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Mxy  (x,  y )  —  (^e^ix  (#1  +  B2  e  2^ix)(y/i  ( Dqq  +  c\  (2  #66  +  ci  #66))  yi  +  ol  ( Dqq  +  c\  (2  #66  +  c\  Hqq 
+  2  Si  Fqq  pi  +  2  ci  £1  #66  Mi)))  + 

0  2m2CC)(^/2  (#66  +  Cl  (2  Fqq  +  Cl  Hqq))  p2  F  OL  (Dqq  +  Cl  (2  Fqq  +  Cl  #66 
+  2  $2  #66  Pi  F  2  Cl  #66  A^)))  + 

(#!  +  #2  e  2>X3X)(^Q  (Dqq  +  Cl  (2  #66  +  Cl  Hqq))  p%  +  QL  (Dqq  +  Cl  (2  #66  +  Cl  #66 
+  2  $3  #66  /i3  +  2  Cl  £3  #66  ^2)))  + 

#4*  (#1  +  #2  e  /i4:E)(74  (#66  +  Cl  (2  #66  +  Cl  Hqq))  p^  F  OL  (Dqq  F  Cl  (2  #66  +  Cl  #66 
+  2(^4  F66  Ml  +  2 Cl  (i.<  B66  Mi))))  cos  (ay)  =  .Mxj,  COS  (ay) 


P;(.x:  (*, y)  =  (eMl  *  (-Bi  +  B2 e“2 x)  (a2  ci  <5i  B12  +  a 71  (Fu  +  ci  B12)  -  Mi  (fif  +  ci  if,,  (1  +  <5i  mi)))+ 

e'"2  *  (_Bl  +  b2  e-^2  (a2  Ci  ^  Hi2  +  a72  (Fl2  +  Cl  ffl2)  -  M2  (Fit  ••  C,  //,,  (1  •  <52  M2)))+ 

eM3  x  (-Bi  +  B2  e_2M3  x)  (a2  ci  <53  H12  +  a  73  (F12  +  ci  Bi2)  -  Ms  (Fu  +  ci  Bn  (1  +  <53  ms)))+ 

eM4  X  (_Bi  +  e-2M4  (q2  Ci  ^  Hi2  +  a74  (Fi2  +  Ci  ffi2)  _  M4  (Fll  +  Cl  Hll  (1  +  <54  M4)))) 


sin  (a  y)  =  Vxx  sin  (a  y) 


(25) 


At  this  point,  zero  boundary  conditions  are  generally  used  to  eliminate  the  constants 
when  using  the  classical  method  which  establishes  the  stability  equation  for  a  single 
individual  plate.  By  contrast,  the  development  of  the  dynamic  stiffness  matrix  entails 
imposition  of  general  boundary  conditions  in  algebraic  form  and  widens  the  possibility 
of  the  analysis  of  multi-plate  systems.  In  order  to  develop  the  dynamic  stiffness  matrix, 
the  following  boundary  conditions  are  applied  next. 


%  —  0  •  fAf m  —  Wmi  5  x m  —  ^xi  5  —  ^yi  1  bbm)£C  —  Wml^x 

X  =  b  :  Wm  =  m2  5  ^Xm  —  <b;i;2  5  ^JJ‘2  5  Wm,x  =  bbm2)X 

(26) 

=  0  •  2a:  =  2a:i  >  ^[xx  ~  J^lxx  1  5  -M-xy  ~  -M-xyi  •>  ^Pxx  ~  ^Pxx± 

X  —  1)  :  Qa:  =  2a?2  )  =  ^PXX2  5  -M-xy  ~  AdXy2  •>  'Pxx  =  ^Pxx2 

By  substituting  Eq.  (26)  into  Eq.(69),  the  following  matrix  relations  for  the  displace¬ 
ments  are  obtained: 


"  Wi  " 

r  5i  -*i 

^2 

-82 

^3 

—  ^3 

64 

-54 

r  Si 

i  i 

1 

1 

1 

1 

1 

1 

s2 

i 

71  -71 

72 

-72 

73 

-73 

74 

-74 

s3 

Wl,« 

fl  ~fi 

/2 

-/2 

fs 

-/3 

/4 

-/4 

ba 

w2 

= 

— 

5i  eb  ^°l  —8\e~b^ol 

eb^o2 

-<52  e~b^o2 

S3  eb^o3 

-<53  e~b^o3 

54eb^o4 

-54e“b^o4 

b5 

$*2 

e^ol  _e-^Bol 

ebB'o2 

_e-f>^o2 

eb  Bo3 

-e~b  ^o3 

eb  A^o4 

_e-b  ^o4 

b6 

71  eb,io1  “71 

^2  ^o2 

—  72  e  —  ^  ^o2 

73 

—  73  e  —  ^  ^°3 

74  Mo4 

—74  e-b 

b7 

w2x 

.  Zie^ol 

f2eblxo2 

~f2e-b^o2 

/3  eb  ^o3 

— f3e-b^o3 

/4e^o4 

-/4e-b^o4  _ 

-  B8 

(27) 

where 

ft 

=  Bi  i 

withi  —  1. 

,2,3,4 

Equations 

(82)  and  (27)  can  be  written 

as 

<5  =  AC  (28) 
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By  applying  the  same  procedure  for  forces  and  moments,  i.e.  substituting  Eq.  (26)  into 
Eq.(25)  the  following  matrix  relations  are  obtained: 


QX1 

MXX1 

VXX1 

Qx2 

mXX2 

j/2 

•MXx2 


where 


Qi 

Qi 

Q2 

Q2 

Qs 

Qs 

Qa 

Q4  1 

-  B1 

Ti 

-Ti 

r2 

-r2 

r3 

-r3 

r4 

-r4 

b2 

-Zi 

-Zi 

-x2 

-x2 

-r3 

-i3 

-x4 

-x4 

b3 

£1 

~£l 

£-2 

—  £ 2 

£ 3 

-£3 

Vl 

-£a 

Ba 

Ql  eb»o  1 

-Ql  e~b^o  1 

Q2  eb^o2 

—  Q2  e~b»o2 

Q3eb^o3 

—  Q3  e~b»o3 

Q4  eb/xo4 

-QAe~b^oA 

b5 

-rieb^oi 

T\e~b  M°1 

-T2eb^o2 

T2e~b»o2 

-T3eb^o3 

T3e~b^o3 

-7^eb^o4 

TAe~b^oA 

b6 

X1eb»ol 

X\e~b ^ol 

X2eb^o2 

Z2e~b^o2 

X3eb»o3 

X3e~b»o3 

X4eb^o4 

X4e~b^oA 

B  7 

—  C±eb  ^°1 

C\e~b  ^ol 

-C2eb^o2 

C2e~b  ^o2 

-C3eb^o3 

C3e~b^o3 

-£4eb/i°4 

C4e~b^oA 

-  B  8 

(29) 


Qi  —  ^55(1  “1“  $ihoi)  2c2(-D55  ~\~  D§§8ifi0i)  ^2(^55  +  8iF^fl0f)  C\  (oL^i  (-Pi 2  F  2 Fqq  F  C\H\2 

+  2ci  H^poi  —  ^{Fn  F  c\H\\  +  ciSiHufioi)  F  o?(2  Fqq  F  2ciHqq  F  c\8iH\2[i0%  F  4ci^i766Moi)) 


Ti  —  o?  c\8i(F\2  +  C1H12)  F  oiyi(Di2  F  ci(2Ei2  +  C1H12))  —  /i oi(Du  +  ci(2En  +  c\H\\ 

+  SiFllhoi  +  CiSiHufloi)) 


F{  —  yi  (Dqq  +  C\(2jFqq  +  CiHqq^ /J.Qi  &(Dqq  F  C\(2jFqq  C\Hqq  F  2SiFQQfX0i  F  2ci(52-f/^66/^oz)) 


A  =  c\[o? c\8{H\2  +  «7i(Fi2  +  C1H12)  -  fJLoi(Fn  F  ciiJn (1  +  SifjLoi)))  with  i  =  1,  2,  3,4 


(30) 


Equation  (29)  can  be  written  as 


F  =  RC 


By  eliminating  the  constants  vector  C  form  Eqs.  (28)  and  (31)  the  dynamic  stiffness 
matrix  is  formulated  as  follows: 

K  =  RA1 

or  more  explicitly 


(31) 
fness 

(32) 


K  = 


Dqq  °qm 
Srnrn 


Sqt 

Sqh 

fqq 

fqrn 

fqt 

fqh 

S'rnt 

Smh 

fqrn 

frnrri 

fmt 

fmh 

Sft 

Sth 

fqt 

fmt 

ftt 

fth 

Shh 

fqh 

fmh 

fth 

fhh 

Sym 

sqq 

Sqm 

Srnrn 

Sqt 

Srnt 

Stt 

- Sqh 

Smh 

- Sth 

Shh 

(33) 
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Finally  the  dynamic  stiffness  matrix  related  to  the  force  and  displacement  vectors  can 
be  written  as  follows: 


Qcc  1 

Sqq  Sqm 

MlXX 1 

Smm 

Mtxyi 

VXX1 

Qx  2 

A/txx  2 

AiXy2 

Vxx2 

Sqt 

Sqh 

fqq 

fqm 

Smt 

Smh 

—  fqrn 

fmm 

Sit 

Sth 

fqt 

fmt 

Sym 

Shh 

~fqh 

Sqq 

fmh 

Sqm 

Smm 

fqt  fqh 

"  Wi  " 

fmt  fmh 

ftt  fth 

fth  fhh 

Wi'X 

Sqt  Sqh 

w2 

Smt  Smh 

Stt  Sth 

®y2 

Shh  _ 

.  ^2,  . 

(34) 


which  in  compact  matrix  form: 

F  =  KD 


(35) 


The  above  dynamic  stiffness  matrix  will  now  be  used  in  conjunction  with  the  Wittrick- 
Williams  algorithm  [13]  to  analyze  composite  simple  and  stiffened  plates  for  their  buck¬ 
ling  behavior  based  on  HSDT.  Explicit  expressions  for  each  element  of  the  DS  matrix 
were  obtained  via  symbolic  computation,  but  they  are  far  too  extensive  and  voluminous 
to  report  here.  The  correctness  of  these  expressions  was  further  checked  by  implementing 
them  in  a  MATLAB  program  and  carrying  out  a  wide  rage  of  numerical  simulations. 


2.4  Assembly  procedure,  boundary  conditions  and  similarities  with 
FEM 

Once  the  DS  matrix  of  a  laminate  element  has  been  developed,  it  can  be  rotated  and/or 
offset  if  required  and  thus  can  be  assembled  to  form  the  global  DS  matrix  of  the  final 
structure.  The  assembly  procedure  is  schematically  shown  in  Fig.  2  which  is  similar  to 
that  of  FEM.  Although  like  the  FEM,  a  mesh  is  required  in  the  DSM,  it  should  be  noted 
that  unlike  the  former,  the  latter  is  not  mesh  dependent  in  the  sense  that  additional 
elements  are  required  only  when  there  is  a  change  in  the  geometry  of  the  structure. 
A  single  DS  laminate  element  is  enough  to  compute  any  number  of  its  buckling  loads 
to  any  desired  accuracy,  which,  of  course,  is  impossible  in  the  FEM.  However,  for  the 
type  of  structures  under  consideration  DS  plate  elements  do  not  have  point  nodes,  but 
have  line  nodes  instead.  Also  no  change  in  the  geometry  along  longitudinal  direction 
is  admitted.  This  assumption  is  in  addition  to  the  assumed  simple  support  boundary 
conditions  on  two  opposite  sides.  The  other  two  sides  of  the  plate  can  have  any  boundary 
conditions.  The  application  of  the  boundary  conditions  of  the  global  dynamic  stiffness 
matrix  involves  the  use  of  the  so-called  penalty  method.  This  consists  of  adding  a  large 
stiffness  to  the  appropriate  position  on  the  leading  diagonal  term  which  corresponds  to 
the  degree  of  freedom  of  the  node  that  needs  to  be  suppressed.  It  is  thus  possible  to  apply 
free,  simple  support  and  clamped  boundary  conditions  on  the  structure  by  penalizing 
the  appropriate  degrees  of  freedom.  Note  that  in  accordance  with  the  notation  and  sign 
convention  used  in  Fig.  2  for  simple  support  boundary  condition  V,  W  and  <f>y  have  to 
be  penalized  whereas  for  clamped  boundary  condition  [/,  V,  W,  W,  x  have  to 

be  penalized.  Clearly  for  the  free-edge  boundary  condition  no  penalty  will  be  applied. 
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Because  of  the  similarities  between  DSM  and  FEM,  DS  elements  can  be  implemented  in 
FEM  codes  and  thus  the  accuracy  of  results  can  be  increased  substantially. 


Dynamic  Stiffness  Matrices 


Figure  2:  Direct  assembly  of  dynamic  stiffness  elements. 


2.5  Application  of  the  Wittrick-Williams  Algorithm 

In  order  to  compute  the  critical  buckling  loads  of  a  structure  by  using  the  DSM,  an 
efficient  way  to  solve  the  eigen-problem  is  to  apply  the  Wittrick  and  Williams  algorithm 
[13]  which  has  featured  in  literally  hundreds  of  papers.  For  the  sake  of  completeness  the 
procedure  is  briefly  summarized  as  follows. 

First  the  global  dynamic  stiffness  matrix  of  the  final  structure  K*  is  computed  for  an 
arbitrarily  chosen  trial  critical  buckling  load  A*.  Next,  by  applying  the  usual  form  of 
Gauss  elimination  the  global  stiffness  matrix,  is  transformed  into  its  upper  triangular 
A*  form.  The  number  of  negative  terms  on  the  leading  diagonal  of  K*  is  now  defined 
as  the  sign  count  s (K*)  which  forms  the  fundamental  basis  of  the  algorithm.  In  its 
simplest  form,  the  algorithm  states  that  j,  the  number  of  critical  buckling  loads  (A)  of 
a  structures  that  he  below  an  arbitrarily  chosen  trial  buckling  load  (A*)  is  given  by: 

j  =  j0  +  s(K*)  (36) 

where  jo  is  the  number  of  critical  buckling  loads  of  all  single  strip  elements  within  the 
structure  which  are  still  lower  than  the  trial  buckling  load  (A*)  when  their  opposite 
sides  are  fully  clamped.  It  is  necessary  to  account  for  this  clamped-clamped  critical 
buckling  loads  because  exact  buckling  analysis  using  DSM  allows  an  infinity  number  of 
critical  buckling  loads  to  be  accounted  for  when  all  the  nodes  of  the  structures  are  fully 
clamped,  (i.e.  in  the  overall  formulation  K  5  =  0,  these  critical  buckling  loads  correspond 
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to  S  =  0  modes.)  Thus  jo  is  an  integral  part  of  the  algorithm  and  is  not  a  peripheral 
issue.  However,  jo  is  usually  zero  and  the  dominant  term  of  the  algorithm  is  the  sign- 
count  s(K*)  of  Eq.  (36).  One  way  of  avoiding  the  computation  of  troublesome  jo  is 
to  split  the  structure  into  sufficient  number  of  elements  so  that  the  clamped-clamped 
buckling  loads  of  an  individual  element  in  the  structure  are  never  exceeded.  Once  s(K *) 
and  jo  of  Eq.  (36)  are  known,  any  suitable  method,  for  example,  bi-section  technique 
can  be  devised  to  bracket  any  critical  buckling  load  within  any  desired  accuracy.  The 
mode  shapes  are  routinely  computed  by  using  the  standard  eigen-solution  procedure  in 
which  the  global  dynamic  stiffness  matrix  is  computed  at  the  critical  buckling  load  and 
the  force  vector  is  set  to  zero  whilst  deleting  one  row  of  the  DS  matrix  and  giving  one  of 
the  nodal  displacement  component  an  arbitrarily  chosen  value  and  determining  the  rest 
of  the  displacements  in  terms  of  the  chosen  one. 


3  Results  ans  Discussion 

Results  are  shown  accounting  for  the  in-plane  load  conditions  in  Fig.  3  and  the  boundary 
conditons  in  Fig.  4.  A  preliminary  validation  of  the  critical  buckling  load  analysis 


Figure  3:  Laminated  composite  plate  under  in-plane  loadings. 

for  moderately  thick  (a/h  =  10)  simply-supported  cross-ply  square  plates  uniaxially 
loaded  in  the  x  direction  is  carried  out  and  the  results  are  shown  in  Table  1  for  different 
orthotropy  ratos.  The  dimensionless  critical  buckling  load,  obtained  using  HSDT  within 
the  framework  of  the  DSM  are  in  excellent  agreement  when  compared  with  the  3D 
elasticity  solution  and  the  results  also  lead  to  the  same  findings  of  the  classical  Levy- 
type  closed  form  solution.  Note  that  for  all  practical  purpose,  it  is  only  the  first  buckling 
load  that  matters.  Therefore  only  the  first  critical  loads  is  presented  in  this  paper.  As 
expected  the  percentage  error,  with  respect  to  the  3D  elasticity  solution,  increases  when 
increasing  the  orthotropic  ratio.  In  Table  2  the  dimensionless  critical  buckling  load  for 
the  same  case  study  of  Table  1  is  computed  but  taking  into  account  the  effects  of  the 
length-to-thickness  and  orthotropy  ratios  and  boundary  conditions  (see  Fig.  4).  At  a 


16 


Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


FA8655-10-1-3084  Report  6  Dynamic  Stiffness  Modelling  of  Plate  and  Shell  Assemblies 


y 


sssc 


'/ 

/ 


X 


Figure  4:  Boundary  conditions. 


Table  1:  Dimensionless  uniaxial  buckling  load  (along  x  direction)  Ncr  =  Ncr  E^h3 ,  for 
simply  supported  cross-ply  square  plates  with  b/h  —  10,  E1/E2  —  open,  G12/E2  = 
G\2t/ E2  =  0.6,  G23/E2  0.5,  v  12  =  v\3  —  0.25. 


Stacking  Sequence 

Models 

E1/E2 

[0o/9079070o] 

3 

20 

40 

3D-Elasticity  [24] 

5.304 

AgD% 

15.019 

A3D  % 

22.881 

A3D% 

Classical  Levy’s  solution 

HSDT 

5.393 

(1.68) 

15.298 

(1.86) 

23.340 

(2.01) 

FSDT 

5.399 

(1.79) 

15.351 

(2.21) 

23.453 

(2.50) 

CLPT 

5.754 

(8.48) 

19.712 

(31.2) 

36.160 

(58.0) 

DSM 

HSDT 

5.393 

(1.68) 

15.298 

(1.86) 

23.340 

(2.01) 

f  A3D%  =  77^-  x  100. 

!  CJQTV 
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fixed  length-to-thickness  ratio,  the  dimensionless  critical  buckling  load  increases  when 
increasing  the  orthotropic  ratio  for  all  the  considered  boundary  conditions.  A  similar 
behavior  can  be  observed  when  varying  the  length-to-thickness  ratio  but  by  fixing  the 
orthotropic  ratio.  Understandably,  the  largest  dimensionless  critical  buckling  load  is 


Table  2:  Dimensionless  uniaxial  buckling  load  (along  x  direction)  Ncr  =  Ncr  /?3 , 
for  simply  supported  cross-ply  square  plates,  stacking  sequence  [0°/90o/90o/0°]  and 
E1/E2  —  open,  G12/E2  —  G13/E2  —  0.6,  G23/E2  —  0.5,  v  12  —  v\z  —  0.25. 


Uniaxial  Compression  (UA-C) 


b/h 


Ei/E2 

BCs 

2 

5 

10 

25 

50 

100 

3 

S-S-S-S 

1.9557 

4.5458 

5.3933 

5.6928 

5.7384 

5.7499 

S-S-S-F 

1.5912 

3.3640 

4.3496 

5.4623 

6.2621 

6.8168 

S-S-S-C 

1.9821 

5.4309 

7.1169 

7.8176 

7.9304 

7.9592 

S-F-S-F 

1.5839 

3.2200 

4.1087 

5.1745 

6.0137 

6.7261 

S-C-S-F 

1.5954 

3.3805 

4.3651 

5.4686 

6.2624 

6.8599 

S-C-S-C 

2.0555 

6.8502 

10.273 

12.012 

12.314 

12.392 

10 

S-S-S-S 

2.2810 

7.1554 

9.9406 

11.209 

11.420 

11.473 

S-S-S-F 

2.0267 

5.1516 

6.9540 

8.6366 

9.7745 

10.536 

S-S-S-C 

2.3426 

8.3895 

14.133 

17.865 

18.587 

18.778 

S-F-S-F 

2.0040 

4.5982 

6.0871 

7.9721 

9.8185 

11.836 

S-C-S-F 

2.0746 

7.1554 

7.0981 

9.0614 

10.724 

12.305 

S-C-S-C 

2.5134 

9.9930 

20.515 

30.068 

32.273 

32.878 

20 

S-S-S-S 

2.5689 

9.4219 

15.298 

18.825 

19.482 

19.654 

S-S-S-F 

2.3144 

6.7586 

9.8649 

12.398 

14.130 

15.513 

S-S-S-C 

2.6643 

10.376 

20.977 

31.075 

33.495 

34.166 

S-F-S-F 

2.2874 

5.8285 

8.1129 

10.893 

13.785 

17.359 

S-C-S-F 

2.4425 

7.0390 

10.416 

13.613 

16.334 

19.282 

S-C-S-C 

2.9945 

11.723 

28.670 

52.336 

59.682 

61.870 

40 

S-S-S-S 

3.0749 

11.997 

23.340 

33.131 

35.347 

35.953 

S-S-S-F 

2.6827 

8.7789 

14.723 

19.374 

22.292 

24.982 

S-S-S-C 

3.1136 

12.301 

29.414 

54.096 

62.190 

64.645 

S-F-S-F 

2.6558 

7.5209 

11.386 

15.744 

20.251 

26.398 

S-C-S-F 

2.9731 

9.4057 

15.987 

21.955 

26.517 

31.952 

S-C-S-C 

3.5992 

13.357 

37.045 

87.668 

110.86 

118.82 

given  by  the  boundary  condition  S-C-S-C  and  the  lowest  by  S-F-S-F.  In  Table  3  the 

18 


Distribution  A:  Approved  for  public  release;  distribution  is  unlimited 


FA8655-10-1-3084  Report  6  Dynamic  Stiffness  Modelling  of  Plate  and  Shell  Assemblies 


results  are  given  for  composite  plates  that  are  uniaxially  loaded  in  the  y  direction,  instead 
of  the  x  direction  for  different  values  of  length-to-t hickness  and  orthotropy  ratios.  The 
dimensionless  critical  buckling  load  is  generally  lower  for  all  the  boundary  conditions 
but  for  for  the  case  with  one  or  two  sides  free,  namely,  S-S-S-F  and  S-F-S-F,  it  decreases 
significantly. 


Table  3:  Dimensionless  uniaxial  buckling  load  (along  y  direction)  Ncr  =  Ncr  h3 , 
for  simply  supported  cross-ply  square  plates,  stacking  sequence  [0°/90o/90o/0°]  and 
E1/E2  —  open,  G12/E2  —  G13/E2  —  0.6,  G23/E2  —  0.5,  v  12  —  v\z  —  0.25. 


Uniaxial  Compression  (UA-C) 


b/h 


Ei/E2 

BCs 

2 

5 

10 

25 

50 

100 

3 

S-S-S-S 

1.7611 

4.5458 

5.3933 

5.6928 

5.7384 

5.7499 

S-S-S-F 

0.9288 

1.4363 

1.5677 

1.6136 

1.6219 

1.6247 

S-S-S-C 

1.8052 

5.0418 

7.0967 

8.0978 

8.2713 

8.3164 

S-F-S-F 

0.6575 

0.9484 

1.0136 

1.0338 

1.0369 

1.0377 

S-C-S-F 

1.0287 

1.7630 

2.0138 

2.1123 

2.1309 

2.1371 

S-C-S-C 

1.8552 

5.3906 

8.0926 

9.6793 

9.9819 

10.062 

10 

S-S-S-S 

1.8421 

6.0365 

9.2387 

10.877 

11.161 

11.234 

S-S-S-F 

1.0769 

1.9529 

2.2269 

2.3238 

2.3402 

2.3452 

S-S-S-C 

1.9124 

6.4505 

10.676 

13.497 

14.069 

14.222 

S-F-S-F 

0.8416 

1.4919 

1.6799 

1.7416 

1.7508 

1.7531 

S-C-S-F 

1.2723 

2.6410 

3.2331 

3.4805 

3.5247 

3.5380 

S-C-S-C 

1.9998 

6.9366 

12.551 

17.665 

18.940 

19.301 

20 

S-S-S-S 

1.9042 

7.2304 

12.728 

16.249 

16.923 

17.100 

S-S-S-F 

1.1977 

2.5793 

3.1263 

3.3311 

3.3647 

3.3741 

S-S-S-C 

2.0019 

7.6623 

14.591 

20.721 

22.214 

22.634 

S-F-S-F 

0.9876 

2.1487 

2.5906 

2.7494 

2.7737 

2.7798 

S-C-S-F 

1.4691 

3.6342 

4.8200 

5.3662 

5.4629 

5.4902 

S-C-S-C 

2.1313 

8.2008 

16.897 

27.735 

31.291 

32.391 

40 

S-S-S-S 

1.9744 

8.5588 

18.034 

26.395 

28.285 

28.801 

S-S-S-F 

1.3172 

3.5410 

4.7817 

5.3152 

5.4037 

5.4272 

S-S-S-C 

2.1210 

8.9972 

20.099 

33.690 

37.996 

39.317 

S-F-S-F 

1.1294 

3.1544 

4.2674 

4.7367 

4.8124 

4.8317 

S-C-S-F 

1.5871 

5.0774 

7.6226 

9.0324 

9.2956 

9.3669 

S-C-S-C 

2.3293 

9.5963 

22.614 

44.455 

54.516 

58.136 
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Part  II 


A  Refined  Dynamic  Stiffness  Element  for  Free  Vibration 
Analysis  of  Cross-Ply  Laminated  Composite  Cylindrical 
and  Spherical  Shallow  Shells 
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Summary 


An  exact  free  vibration  analysis  of  laminated  composite  doubly-curved  shallow  shells 
has  been  carried  out  by  combining  the  dynamic  stiffness  method  (DSM)  and  a  higher 
order  shear  deformation  theory  (HSDT)  for  the  first  time.  In  essence,  the  HSDT  has  been 
exploited  to  develop  first  the  dynamic  stiffness  (DS)  element  matrix  and  then  the  global 
DS  matrix  of  composite  cylindrical  and  spherical  shallow  shell  structures  by  assembling 
the  individual  DS  elements.  As  an  essential  prerequisite,  Hamilton’s  principle  is  used  to 
derive  the  governing  differential  equations  and  the  related  natural  boundary  conditions. 
The  equations  are  solved  symbolically  in  an  exact  sense  and  the  DS  matrix  is  formulated 
by  imposing  the  natural  boundary  conditions  in  algebraic  form.  The  Wittrick-Williams 
algorithm  is  used  as  a  solution  technique  to  compute  the  eigenvalues  of  the  overall  DS 
matrix.  The  effect  of  several  parameters  such  as  boundary  conditions,  orthotropic  ratio, 
length-to-thickness  ratio,  radius-to-length  ratio  and  stacking  sequence  on  the  natural 
frequencies  and  mode  shapes  is  investigated  in  details.  Results  are  compared  with  those 
available  in  the  literature.  Finally  some  concluding  remarks  are  drawn. 
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4  Introduction 

Aerospace  structures  are  generally  made  up  of  thin-walled  cylindrical  or  spherical  shell 
components.  The  application  of  composite  shell  structures  is  justified  because  of  the 
extraordinary  load-carrying  capability.  Furthermore,  their  structural  efficiency  is  char¬ 
acterized  by  a  high  stiffness-to-weight  and  strength-to-weight  ratios.  Due  to  the  harsh 
environment  conditions  in  which  aerospace  structures  work  a  comprehensive  and  accu¬ 
rate  investigation  of  their  dynamic  behavior  is  required.  During  the  last  decades  many 
theories  have  been  developed  to  this  purpose.  In  particular,  most  of  the  classical  theo¬ 
ries  are  mainly  based  on  Kirchhoff-Love’s  hypotheses  [25,26].  It  is  well-known  that  these 
shell  theories  neglecting  the  transverse  shear  and  normal  stresses  lead  to  accurate  results 
when  high  values  of  length-to-thickness  ratio  are  considered  and  when  3D  local  effects 
are  not  present.  When  moderately  thick  or  thick  shells  are  investigated  the  inclusion 
of  transverse  shear  and  normal  stresses  are  mandatory.  However,  on  the  other  hand,  a 
refinement  in  results  accuracy  cannot  be  only  achieved  by  virtue  of  an  enhancement  in 
the  kinematics  description  of  the  displacement  model,  but  combining  it  with  an  adequate 
description  of  the  curvature  terms  h/Ri  with  i  =  cq  /3.  In  this  regards,  some  interesting 
observations  have  been  provided  by  Qatu  [27]  and,  Carrera  [28, 29]  and  recently  Fazzo- 
lari  [5,30].  In  order  to  provide  solutions  of  the  Governing  Differential  Equations  (GDEs) 
of  practical  interest,  considerable  efforts  were  focused  and  devoted  to  their  simplifica¬ 
tion.  As  a  result,  the  GDEs  of  shallow  shells  were  derived.  Most  notably,  it  worths 
highlighting  the  articles  of  Donnell  [31,32],  Mushtari  [33,34]  and  Vlasov  [35,36].  The 
hypothesis  of  thinness  was  discarded  in  the  shell  theories  proposed  by  Fliigge  [37,38], 
Lur’e  [39]  and  Byrne  [40],  where  higher-order  expansion  of  the  reciprocal  of  the  Lame 
parameters  was  considered.  Other  refinements  of  the  developed  shell  theories  were  pro¬ 
posed  by  Sanders  [41].  Additional  effects  in  the  development  of  shell  theories  were  taken 
into  account  by  Whitney  and  Sun  [42],  Librescu  [43],  Gulati  and  Essemberg  [44],  Zucas 
and  Vinson  [45]  and  Ambartsumian  [46-53]  amongst  others.  Additional  references  can 
be  found  in  Naghdi  [54],  Ambartsumian  [55]  and  Bert  [56-58].  Reddy  [59]  proposed  a 
generalization  of  Sander’s  theory  to  anisotropic  doubly-curved  shells.  The  application 
of  layer-wise  theories  for  shell  structures  can  be  found  in  the  papers  presented  by  Hsu 
and  Wang  [60],  Cheung  and  Wu  [61],  Barbero  et  al.  [62]  and  Carrera  [28,29,63].  Re¬ 
views  on  finite  element  shell  formulations  have  been  given  by  Denis  and  Palazzotto  [64] 
and  Di  and  Ramm  [65].  An  exhaustive  review  on  classical  theories  can  be  found  in 
Librescu  [43].  As  regards  the  use  of  approximation  methods,  Qatu  and  Asadi  [66]  ad¬ 
dressed  the  vibration  analysis  of  doubly-curved  shallow  shells  with  arbitrary  boundary 
conditions  by  using  the  Ritz  method  with  algebraic  polynomial  displacement  functions. 
Asadi  et  al.  [67]  employed  a  3D  and  several  shear  deformation  theories  in  order  to 
carry  out  static  and  vibration  analysis  of  thick  deep  laminated  cylindrical  shells.  Fer¬ 
reira  et  al.  [68]  used  a  wavelet  collocation  method  for  the  analysis  of  laminated  shells. 
The  same  author  [69]  combined  a  sinusoidal  shear  deformation  theory  with  the  radial 
basis  functions  collocation  method  to  deal  with  static  and  vibration  analyses  of  lami¬ 
nated  composite  shells.  Tornabene  et  al.  [70,  71]  studied  the  free  vibration  behavior  of 
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doubly-curved  anisotropic  laminated  composite  shells  and  revolution  panels  by  means 
of  the  Generalized  Differential  Quadrature  (GDQ).  Fazzolari  et  al  [5,30,72-76]  pro¬ 
posed  two  advanced  hierarchical  trigonometric  Ritz  formulations  based  on  the  Principle 
of  Virtual  Displacements  (PVD)  and  Reissner’s  Mixed  Variational  Theorem  (RMVT)  to 
cope  with  the  free  vibrations  of  anisotropic  laminated  composite  plates  and  cylindrical, 
spherical  and  hyperbolic  paraboloidal  shells.  The  dynamic  stiffness  method  (DSM)  has 
not  been  extensively  researched  in  the  analysis  of  laminated  composite  shells,  only  few 
works  have  been  devoted  to  the  application  of  DSM  to  shell  structures  and  very  few  to 
laminated  composite  cylindrical  and  spherical  shells.  In  particular,  Casimir  et.  al  [77] 
developed  the  dynamic  stiffness  matrix  of  tubular  shells  with  free  edge  boundary  con¬ 
ditions,  including  rotatory  inertia,  transverse  shear  deformation  and  non-axisymmetric 
loadings,  but  the  analysis  was  restricted  to  homogeneous  and  isotropic  material.  Thinh 
and  Nguyen  [78]  proposed  the  dynamic  stiffness  matrix  of  continuous  element  for  vibra¬ 
tion  of  thick  cross-ply  laminated  composite  cylindrical  shells.  Khadimallah  et.  al  [79] 
dealt  with  the  dynamic  response  to  harmonic  loads  of  an  axisymmetric  shell  by  using 
a  dynamic  stiffness  formulations.  The  dynamic  stiffness  matrix  of  toroidal  shells  was 
derived  by  Leung  and  Kwok  [80].  El-Kaabazi  and  Kennedy  [81]  proposed  dynamic  stiff¬ 
ness  equations  for  variable  thickness  cylindrical  shells  under  the  assumptions  of  Donnell, 
Timoshenko  and  Fliigge  theories.  Langley  [82]  developed  a  dynamic  stiffness  technique 
for  the  vibration  analysis  of  stiffened  shell  structures.  Chronopoulos  et.  al  [83]  by  using 
a  dynamic  stiffness  approach  studied  the  response  of  layered  shells.  Tounsi  et.  al  [84]  de¬ 
rived  the  dynamic  stiffness  matrix  of  circular  rings.  In  the  present  article  a  new  dynamic 
stiffness  element  has  been  developed  in  order  to  cope  with  the  free  vibration  analysis  of 
laminated  composite  doubly-curved  shallow  shells.  The  proposed  formulation  is  based 
on  the  close-form  Levy-type  solution  and  the  application  of  the  Wittrick-Williams  al¬ 
gorithm  [85].  The  latter  is  used  as  a  solution  technique  to  compute  the  eigenvalues 
of  the  overall  DS  matrix.  Several  symmetric  cross-ply  lamination  schemes  have  been 
investigated.  The  effect  of  different  parameters  such  as  stacking  sequence,  boundary 
conditions,  orthotropic  ratio,  length-to-thickness  and  radius-to-length  ratios  on  the  di¬ 
mensionless  frequencies  parameters  and  mode  shapes  have  been  studied.  Finally  some 
conclusions  have  been  drawn  from  the  findings  of  the  research. 

5  Theoretical  Formulation 

In  the  derivation  that  follows,  the  hypotheses  of  straightness  and  normality  of  a  trans¬ 
verse  normal  after  deformation  are  assumed  to  be  no  longer  valid  for  the  displacement 
field  which  is  now  considered  to  be  a  cubic  function  in  the  thickness  coordinate  and 
hence  the  use  of  higher  order  shear  deformation  theory  (HSDT).  The  composite  circular 
cylindrical  shell  is  assumed  to  be  composed  of  Ni  layers  so  that  the  theory  is  sufficiently 
general.  The  geometry  of  the  shell  and  the  coordinate  system  are  shown  in  Fig.  5. 
The  integer  k  is  used  as  a  superscript  denoting  the  layer  number  where  the  numbering 
starts  from  the  bottom.  After  imposing  the  transverse  shear  stress  homogeneous  con¬ 
ditions  [19,86]  at  the  top/bottom  surface  of  the  shell,  the  displacements  field  is  in  the 

23 


Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


FA8655-10-1-3084  Report  6  Dynamic  Stiffness  Modelling  of  Plate  and  Shell  Assemblies 


usual  form  given  below: 


u  ( a ,  P,  z,  t) 


v(a,P,z,t) 
w  ( a ,  P,  z,  t) 


(l  +  y'j  Mo  (a,  P  ,t)+z  4>a  (a,  P,  t )  -  z 3  ^ 

(1  +  ifs)  V°  P,t')  +  Z  P’ ““  z*  3^2 
wo  ( a,P,t ) 


(<£c«  (a,P,t)  + 
[<t>(3  (a,P,t)  + 


1  dwo{a,p,t)\ 
Aa  da  ) 

1  dwo(a,P,t)\ 

A p  dp  J 


(37) 


where  u,  v,  w  are  general  displacements  within  the  shell  in  the  <a,  /?,  and  z  directions, 
respectively,  whereas  uq,  vo ,  u;o  are  the  corresponding  displacements  of  the  reference 
surface  (mid-plane  FI),  0a,  are  the  rotations  of  the  normals  to  the  cq  /?,  respectively, 
Aa,  Ap  are  the  surface  metrics  and  i?a  and  Rp  are  the  radii  of  curvature  in  the  a  and  (3 
directions.  The  geometry  of  the  shell  is  completely  described  by  the  parameters  given  in 
Fig.  5.  Most  notably,  r  (a,  ft)  is  the  position  vector  of  a  point  on  the  middle  surface  Ft  of 


the  shell,  R  (a,  /?,  z)  is  the  position  vector  of  a  generic  point  within  the  volume  occupied 
by  the  shell.  At  each  point  P  of  the  middle  surface  n  (a,  (3)  indicates  the  unit  normal 
vector.  Subscripts  and  superscripts  k  refers  that  quantity  at  layer  level.  The  square  of 
line  elements  in  orthogonal  curvilinear  coordinates  is  therefore  defined  as: 

dsl  =  drk  •  drk  =(H*y  dal  +  (vf)  "  d/3%  +  (tf*)  2  dz\  (38) 

the  area  of  an  infinitesimal  rectangle  on  Fl^  as: 

dnk  =  H*Hkpdakdf3k  (39) 
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and  the  infinitesimal  volume  as: 

dVk 

where 


TTk  _  Ak 


1  +  ^) 

RkJ 


HkHk0Hkdakd3kdzk 


Hi  =  4 


m  =  i 


(40) 


(41) 


are  referred  to  as  Lame  parameters  and  A k  and  A k  are  the  first  fundamental  magnitude 
of  the  first  fundamental  form  dfik-  Attention  is  herein  focused  on  shells  with  a  constant 
curvature,  i.e.,  doubly-curved  shells  (cylindrical,  spherical,  toroidal  geometries)  for  which 
=  Ak3  =  1.  In  the  case  of  a  swallow  shell  the  approximation  (l  + 


Rcy 


1  and 


(1  +  ^) 


1  is  generally  valid  and  leads  to  reasonably  reliable  results.  The  strain- 
displacement  relationships  referred  to  an  orthogonal  curvilinear  coordinate  system  give 
rise  to  the  following  deformation  field  [19,68]: 


£aa  ~  £aa  +  z  (C  +  ^  kao)  >  £PP  ~  ^fBp  +  z  ^  ^/?)  > 

7 a/3  =  7 a0  +  z  (^a/3  +  ^  ^4/3  )  >  7ca  =  7 az  +  ^  4L>  Ipz  =  7/L  +  ^  ^pz 


(42) 


where 


o  9u0 
£aa  ~  da  ' 
o  dw0 

7«=  aV 

0  _  9ffa 
“  da  ’ 


£PP  - 


w 

Ra 


d/3  > 


du0 

93 


w 

Re 


1% 


du0 

93 


9y_  o 
3a 


o  _  9wo  ,  , 

l0z~w  ^ 


uO  =  9<j>a 

af}  d3  +  3a 


,Z  3  9(f)  Q,  92w  o 

= C1 1  + 


,2  _  [  "Ti; 

KB8  ~  °1 


d(f>B  ,  d2w0 
3/3  +  3/32 


k2  =  (93i  93a  9yg\  ! 

pa  1  \  da  +  d3  +  dad3 )  ’ 


(43) 


■c3(^+S’  fc»=c2(^+i) 


The  stress-stain  relationship  of  the  k-layer  is 


ak  =  Ckek 


(44) 


where  C  is  the  plane-stress  constitutive  matrix  (see  [19])  and  ak  and  £  '  are  the  stress 
and  strain  vectors,  respectively,  and  defined  as 


<7*  = 


£fc  = 


{ 


aa 

k 


a 


PP  aap 


Tk 

az 


aa 


’pZ  } 

£pp  7 Ip  ilz  }T 


(45) 


where  T  denotes  a  transpose. 
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5.1  Governing  differential  equations 

The  governing  differential  equations  (GDEs)  are  derived  by  using  Hamilton’s  Principle. 
The  variational  statement  can  be  written  as 


JX  rl 2 

8Ckdt  =  0  (46) 

k= 1  tl 


where  Ck  is  the  Lagrangian  for  the  kth  layer  of  the  composite  shell.  The  first  variation 
can  be  expressed  as 

8Ck  =  8Tk  -  8Uk  (47) 

where  8Uk  is  the  virtual  strain  energy,  8Tk  is  the  virtual  kinetic  energy,  and  assume  the 
following  form 


SUk  = 


r  rzk+1 

/  /  UekT  ak\  dnk  dz ; 

J  Qk  J  zk  '  ' 


(48) 


8V 


P  PZk+1 

=  ( 
Jzk  V 


pk  8rjT  t)  )  dflk  dz 


where  stresses  (<rfc)  and  strains  (efc)  are  defined  in  Eq.  (45)  and  r]  is  the  displacement 
vector  given  by 

rj  =  {  u  v  w  }T  (49) 

pk  denotes  mass  density  while  an  over  dot  denotes  differentiation  with  respect  to  time. 
The  symbol  8  represents  the  variational  operator.  Imposing  the  condition  in  Eq.  (46), 
the  equations  of  motion  are  derived 


8u0  : 
8v o  : 
8w0  : 


dNaa  dNap  _  d2u0  -  d2pa  d  ( d2w0\ 

da  d/3  0  dt 2  1  dt 2  3  da  \  dt2  J; 

dNap  dNpp  d2v o  -  d2<j>p  d  fd2w 0\ 

da  dp  0  dt- 2  1  dt 2  3  dp  \  dt2  )' 


dQ 


da 


a  a  .  dQ 

a  3  . 

+  ^  +C2 


Nn 


dp 


dKpu 

da 


+ 


dK, 


00 


dp 


Cl 


d2Pa 


da2 


+ 


d2P, 


00 


Ra  R/3 

y  ( <yty\ 

dp  \  dt 2  J  J 


=  -ci  h 


d  ^  ®  ( ^2y[ 


dp2 

r  d 


+  2- 


02P, 


a  (3 


da  (8 


-  ci  h  vr- 


,  T  d2w°  2  t 

+  h^P~Clh 


da  V  dt2  )  dp  V  dt2 

d2  ( d2wo\  d2  ( d2wo\ 
da2  \  dt 2  )  dp2  \  dt2  J 


d2p0 


Ida  \  dt 2 


+ 


(50) 
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84>0 


<% 


9Ma  a  ^  dMa  p 


I* 


da 

d2cf>0 


dp 

W  9  fd2W0 


C2  P-ot  q,  +  C\ 


2  dt 2  +  Cl  ^  da  I  dt 2 


dMpa  dMpp 


da 

d2c 


+ 


dp 

r  d  fd2W0 


T *  w  Y 'P  ,  f  ^ 

2  ~Wr  +  Cll4dp  gdP 


The  natural  boundary  conditions 


are 


dPqq  dPq£ 
da  d/3 


Qpp-C2Kpp+Cl  I  ■ 


dPnR  dP , 


+ 


PP 


dp 


?  5V) 


t  i 

h  dt 2  + 


Su0  : 
5v o  : 


8w0  : 

8pa  : 
<%  : 
8w0, a  : 
8w0,p  : 


Nn  —  na  Na  a  +  Tip  Np  p  +  2  nQ  Tip  Na  p 

Nns  =  (P^i3  /?  Not  a)  Hf)  P  Na  p  (jTa  Tl p ) 


q  —  jy  +  N  ®W° 

Vn  —  iVn  o  i  iVns 


dn 


ds 


ci 


dPn  dPns 


dn 


ds 


+  Qn  +  c2  Kn 


Mn  —  Mn  +  Ci  Pn 
M-ns  Afns  T  Ci  Pns 

Pn  —  C\  Pji 


Cl  Pns 


where 

Pn  =  n2a  Paa  +  Ppp  +  2nanp Pap 
Pns  =  ( -P/3/3  -  ^aa)  %  Tip  +  Pa  p  (n2  -  Tip) 

PI 'a  =  Tla  Plot  a  P  Tip  Mp  p  -\-  2  T la  Tip  Ma  p 
Mns  =  ( Mpp  -  Maa)nanp  +  Map  (n2  -  rip ) 
Qn  =  Qaa  Cl a  T  Q (3  /3  Cl (3 
Pn  ~  Pot  a  Clot  T  P/3  (3  Cl (3 


(51) 


(52) 


and  ci  =  —  ,  C2  =  —  p-.  Moreover,  na,  np  are  the  direction  cosines  of  the  unit  normal 

on  the  boundary  of  the  laminate.  The  force  and  moment  resultants  appearing  in  Eqs. 
(50)  and  (51)  are  defined  as  follows: 


Nl  +  l 


(Ni,  Mi,  Pi)  —  /  (of  (l,  Zj  z3)')  dz  with  i  —  aa,  /?/?,  /3a 

k= 1  Jzk  V 

Nl  „zk+ 1 

(Qaa,Kaa)  =  J2  /  (<rL(M2))dz 

fc=i 

AT/  pzk+1 

(Qpp,Kpp)  =  J2  j  (okpzyz2))dz 


(53) 


fc=i 
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(54) 

(55) 


The  inertia  terms  in  Eq.  (50)  are  defined  as 

li  =  li  +  C1  Ii+ 2 ;  I2  =  ^2  +  2  Cl  /4  +  /6 

and 

Afz  nzk~ir^ 

(Io,h,h,h,h,h)  =  T  /  p  (l  z,  z2,  z3,  z4,  z6)  cZz 

fe=i 

The  resultants  in  Eq.  (53)  cab  be  expressed  in  terms  of  strain  components 

{AT}  =  [A]  p0}  +  [5]  p°}  +  [£]  p2} 

{M}  =  [B]  {e°}  +  [D]  {fc°}  +  [F]  { k 2} 

m  =  p]  p°}  +  pi  {fc°}  +  [^]  { k 2} 

{Qi  =  [as]  {70}  +  m  { k 4} 

{K}  =  [DS\  {7°}  +  [Fs]  { k q 

Where  laminate  stiffnesses  are  defined  as 

Ni  oZk+1 

{Aij ,  Bij ,  Dij ,  Eij,Fij ,  Hij)  =  W  /  Cf)  (1  z,  z2,  z3,  z4,  z6)  dz  with  i,j  =  1,  2, 6 

fc=i 

A^z  pzk+1 

C AijiDtj >  *$)  =  E  /  6  z2>  ^4) with  *>  j  = 4’ 5 

7„ _ i  */zfe 

(57) 


(56) 


fc=l  ' 


The  resultants  in  Eq.  (56)  are  collected  as  follows 

=  {-^aa?  ^/3/3 ?  5  =  M Map}  ]  =  -f/3/3?  Pa [3 }  5 

{Q}  =  {Qaz,  Q/3z}  5  {^}  =  {^az, 

(58) 

similarly  for  the  strains  which  are  given  as 

P0}  =  PL>  4P’  7^}  ;  {fc°}  =  {C*,  *$/?,  ;  {k2}  =  {klco  fcj^}  *, 

{7°}  =  {tL,7P};  P1}  =  PL,^} 

(59) 

For  cross-ply  stacking  sequences  an  exact  solution  can  be  sought  in  Navier’s  or  Levy’s 
form.  For  such  laminates  C\§  =  C^q  =  C\ 5  =  0.  Moreover,  if  the  stacking  sequence  is 
symmetric  the  coupling  elastic  coefficients  B \j  and  E{j  (see  Eq.  (57))  are  zero. 
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5.2  Dynamic  Stiffness  Method 

Once  the  equations  of  motion  Eqs.(50)  and  the  natural  boundary  conditions  Eqs.(51) 
are  derived,  the  classical  procedure  to  carry  out  exact  free  vibration  analysis  of  the 
present  type  consists  of  (i)  solving  the  system  of  differential  equations  in  Navier  or  Levy 
type  closed  form  in  an  exact  sense,  (ii)  applying  particular  boundary  conditions  on  the 
edges  and  finally  (iii)  obtaining  the  frequency  equation  by  eliminating  the  integration 
constants  [20-23].  This  method,  although  extremely  useful  for  analyzing  an  individual 
element,  it  lacks  generality  and  cannot  be  easily  applied  to  complex  structures  for  which 
researchers  usually  resort  to  approximate  methods  such  as  the  FEM.  In  this  respect,  the 
dynamic  stiffness  method  has  no  such  limitations  as  it  always  retains  the  exactness  of 
the  solution  even  thought  when  it  is  applied  to  complex  structures.  This  is  because  once 
the  dynamic  stiffness  matrix  of  a  structural  element  is  developed,  it  can  be  offset  and/or 
rotated  and  assembled  in  a  global  DS  matrix  in  the  same  way  as  the  FEM.  This  global 
DS  matrix  contains  implicitly  all  the  exact  natural  frequencies  of  the  structure  which 
can  be  computed  by  using  the  well  established  Wittrick-Williams  algorithm  [85]. 

A  general  procedure  to  develop  the  dynamic  stiffness  matrix  of  a  structural  element  can 
be  summarized  as  follows: 

(i)  Seek  a  closed  form  analytical  solution  of  the  governing  differential  equations  of 
motion  of  the  structural  element  under  consideration  in  free  vibration. 

(ii)  Apply  a  number  of  general  boundary  conditions  that  are  equal  to  twice  the  num¬ 
ber  of  integration  constants;  these  are  usually  nodal  displacements  and  forces  in 
algebraic  forms. 

(iii)  Eliminate  the  integration  constants  by  relating  the  amplitudes  of  the  harmonically 
varying  nodal  forces  to  those  of  the  corresponding  displacements  which  essentially 
generates  the  frequency- dependent  dynamic  stiffness  matrix,  providing  the  force- 
displacement  relationship  of  the  nodal  lines. 


5.3  DS  formulation  and  Levy-type  closed  form  exact  solution 

The  solution  of  the  GDEs  is  sought  in  the  following  form: 


u°(a,/3,£)  =  ^2  Um(ot)  elujt  sin(<9  /?); 

m=  1 
oo 

w°(a,/3,t)  =  *^2  fEm(a)  elujt  sin(0/3); 

m= 1 

oo 

<j>p{a,f},t)  =  y  pm(a)  elut  cos(6 P) 

m=  1 


v°(a,P,t)  =  y  Vm(a)  eluJt  cos(0 p); 

m= 1 
oo 

<t>a(a,P,t)  =  y  $am(a)  elult  sm(6 (3); 

m=  1 


(60) 


where  uo  is  the  unknown  circular  frequency,  6  —  rn^L  and  m  =  1,2,. . . ,  oo.  This  is  the 
so-called  Levy’s  solution  which  assumes  that  two  of  the  opposite  sides  of  the  shell  are 
simply  supported  (S-S),  i.e.  vq  =  wq  —  (f)a  —  0  at  /?  =  0  and  /3  —  b.  Substituting  Eq. 
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(60)  in  the  GDEs  a  set  of  ordinary  differential  equations  (ODEs)  is  obtained  which  can 
be  written  in  matrix  forms  as  follows: 


£11  £12  £13  £14  £15 

£21  £22  £23  £24  £25 

£31  £32  £33  £34  £35 

£41  £42  £43  £44  £45 

£51  £52  £53  £54  £55 


'  0 
0 

=  0 
0 
0 


where  Cij  with  i,  j  =  1,  •  •  •  ,6  are  differential  operators  following  defined: 
£11  =  ^66  a2  —  An  'Dq  —  h  u2-,  £12  =  (^12  +  ^66)  ol  Va; 

Ca  =  (~'4n  “  Al2  v°' 

£14  =  £15  =  £24  =  £25  =  £41  =  £42  =  £51  =  £52  =  0; 

£21  =  £12;  £22  =  “^22  «2  +  ^66  fd2  +  In  lo 2;  £23  =  a  (  Ai?  — — h  An 


£31  =  £13; 


£32  =  £23; 


£33  —  — +  Io  to2  +  2?2  ^55  +  2  C2  -D55  +  c2  E55  +  2  a2  cf  H12  +  4  a2  c2  Hg 
—  cf  Iq  uj2^J  —  a2  (^44  +  C2  (2  D44  +  C2  F44)  +  cf  (a2  H22  —  h  w2))  —  An 


Rn  Rr- 


£34  =  (—ci  Fn  —  cf  Hu)  V ^  +  ^55  +  2  C2  D55  +  a2  ci  F12  +  cf  F55  +  2a2  ciFee 


+  a2  cf  H12  +  2a2  cf  Hqq  -  ci  h  u2  -  cf  J6  J2  )  Va ; 


£35  —  R*a  C1  ^12  +  2  a  Cl  ^66  +  «  C2  H\2  +  2o;c|  #66)  —  a  ^44  +  C2  (2  F44  +  C2  F44); 

+  a2  ci  (F22  +  ci  H22)  ~  c\  (I4  +  ci  Iq)  cu2^j ; 

£43  =  -£34;  £53  =  -£35;  £54  =  —£45; 

£44  =  —^55  ~  C2  (2  £>55  +  C2  F55)  +  R*a  (£*n  +  2  ci  Fn  +  c2  #n)  —  <a2  ^F66  +  2  ci  F66 
+  cf  +  ^/2  +  2  Cl  Z4  +  c2  u2] 

£45  =  (— <a  D12  —  a  Dm  —  2  a  ci  F12  —  2a  ciFqq  —  a  cf  H12  —  a  cf  Hqq)  ; 

£55  =  —^44  —  C2  (2  F44  +  C2  F44)  —  a2  ( D22  +  Cl  (2  F22  +  Cl  #22))  +  Ra  (^66  +  Cl 
(2  _F|36  +  Cl  Hqq))  +  (/2  +  Cl  (2  I4  +  Cl  /g))  Cj2 

(62) 


where  Va  —  ci  and  C2  have  already  been  given  previously,  and  Aij,  Bij ,  CV/,  F^j, 
F^-,  Fij,  Hij  have  already  been  defined  in  Eq.  (57).  Expanding  the  determinant  of  the 
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matrix  in  Eq.  (61)  the  following  ordinary  differential  equation  is  obtained: 

{cl\  D^  +  &2  D^  T  D^  +  &4  Da  Da  T  D 2  +  6Z7)  4/  =  0  (63) 

where 

^  =  Um  or  Vm  or  Wm  or  <ha  or  (64) 

Using  a  trial  solution  e^A  in  Eq.  (63)  yields  the  following  auxiliary  equation: 

cl\  A12  +  $2  A10  +  U3  A8  +  ^4  A6  +  ^4  A4  ~\~  clq  A2  H-  ^7  =  0  (65) 

the  polynomial  coefficients  aj  with  j  =  1,  •  •  •  ,  7  have  been  derived  in  symbolic  form.  The 

details  are  not  reported  here  for  brevity.  Numerically  some  roots  of  Eq.  (65)  may  turn 
out  to  be  complex.  As  complex  roots  occur  in  conjugate  pairs,  the  associated  coefficients 
will  also  occur  in  conjugate  pairs.  The  solution  of  the  system  of  ordinary  differential 
equations  can  thus  be  written  as: 

12  12  12 

Um  («)  =  J2  M  eAi“;  Vm  (a)  =  ]T  eAi Wm  (a)  =  eAi  Q; 

i=1  i=1  i=1  (66) 

12  12  V  J 

$a(a)  =  J]^eA*“;  ^  (a)  =  eAi“ 

2=1  2=1 

where  3^,  <5^,  ^  with  i  =  1,  •  •  •  ,12  are  integration  constants.  Note  that  as  a 

results  of  Eq.  (61)  the  sets  of  the  constants  are  not  all  independent.  Thus  a  set  of  only 
twelve  independent  constants  can  be  chosen  and  they  can  be  related  to  the  other  sets. 
This  choice  is  completely  arbitrary.  In  the  present  case  the  set,  is  chosen  to  be  a 
base.  By  substituting  Eqs.  (66)  into  (61)  the  following  relationships  can  be  obtained: 

~  'Hi  ^2  =  l[i  Xi  -®2  5  ~  & 2  ^2  =  ^2  X2  5  ^2  =  X2  ^ 2  5  ~  02  ^ 2 

(67) 

where 

A i  (A^2Q:2 Rot  +  Tin  (— A22A2  +  loco2  +  AqqX2)  Rq,  +  A22A.QQQ!2 Rp  +  A\2  (Aqqcx,2 Rcx  +  Iqw2  Rp  +  A66X2Rp)) 
( A22C y2  —  low2)  (Aqqo2  —  Iqw2)  +  ((A22  —  A11A22  +  2Ai2Aqo)  a2  +  (An  +  Aqq)Iow2)  X2  +  AnAooXf 


ol  (— (A12  +  Aqq) X2 (An Rex  +  A±2Rp)  —  {Aq qol2  —  low2  —  A±±X2)  (An-Ra  +  A22Rp)) 
(. Ai2aXi  +  AooaXi)2  —  ( Aqq(x 2  —  Iqw 2  —  An  A?)  (— A22a2  +  Iqw2  +  AqqX2) 


(68) 


31 


Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


FA8655-10-1-3084  Report  6  Dynamic  Stiffness  Modelling  of  Plate  and  Shell  Assemblies 


Xi  =  —  (  —  O'2  (£>12  +  £>66  +  Ci(2£i2  +  2£66  +  Cl  (#12  +  Hqq)))2  a2  +  ^  —  A55  —  C2(2£>55  +  C2£55)  —  O?  (£>66  +  Cl  (2£66 
+  c\Hqq))  +  (£2  +  ci (2/4  +  cile))^2  +  (£>11  +  ci(2£\i  +  ci££n))A2^)  (A44  +  02(21)44  +  C2-P44)  +  ck2(£>22 
+  Ci(2£22  +  C1H22))  —  (h  +  Cl  (2/4  +  Ci/6))cj2  —  (£>66  +  Ci(2£66  +  CiHqq))X2^  /  ^  —  Xi  (A44  +  C2(2£>44 
+  C2F44)  +  O?  (£>22  +  Ci(2£22  +  C1H22))  —  (£2  +  Cl  (2/4  +  C\Iq))u2  —  (£>66  +  Ci(2£66  +  CiHqq))X2^  (A55 
+  2c2£>55  +  ^£55  +  Cl  (a2 (£12  +  2£66  +  Ci(ifi2  +  2Hqq))  —  (£4  +  c\Iq)w2  —  (£11  +  cii£n)A2)  ^  +  o? (£>12 
+  £>66  +  ci(2£\2  +  2£66  +  ci(ifi2  +  HQe)))\i  [A44  +  2c2£>44  +  c|£i4  +  ci  (a2(£22  +  C1H22)  —  (£4  +  ci£e)u;2 

—  (£12  +  2£66  +  Cl  (££12  +  2Hqq))X2^^ 

^  =  _a(£>i2  +  066  +  ci(2Fi2  +  2 F66  +  Cl{H12  +  H66)))Xi  ^  +  2C2l>55  +  “  °66  +  +  2“  ClFe6 

+  o?  c2Hqq  —  I2UJ2  —  2c\l4Uj2  —  c2Iqll>2  —  £>11  A2  —  2ci£nA2  —  c2££n  X2  —  (a2(£>i2  +  £>66  +  ci(2£i2  +  2£66 

+  Ci(££i2  +  Hqq)))2X^  (  —  A55  —  C2 (2£>55  +  C2£55)  +  Cl  (  —  a2(£i2  +  2£66  +  ci(££i2  +  2Hqq))  +  (£4  +  ci£e)w2  +  (£11 

+  ci££n)A2^^/^  —  A i  (A44  +  C2 (2£>44  +  C2F44)  +  a2(£>2 2  +  ci(2£22  +  C1H22))  —  (£2  +  ci(2£4  +  ciIq))u2 

—  ( £>66  +  Cl  (2Fqq  +  CiHqq))X2^J  (A55  +  2c2£>55  +  ^£55  +  Cl  [o?  (F12  +  2£66  +  Cl  (££12  +  2Hqq))  —  (£4  +  ci£6)cj2 

—  (£11  +  Ci££n)A2^  +  a2(£>12  +  £>66  +  Ci(2£i2  +  2£66  +  Ci(££i2  +  Hqq))) Xi  (A44  +  2c2£>44  +  C2F44  +  Cl  ^a2(£22 

+  Cl ££22)  —  (£4  +  ciIq)uj2  —  (£12  +  2£66  +  Cl  (££12  +  2Hqq))X2'J'J'J  —  (a i  ^55  +  C2{2Dqq  +  C2F55)  +  a2 (£>66  +  ci(2£66 
+  ci ££66))  —  (£2  +  ci(2£4  +  ciIq))uj2  —  (Du  +  ci(2£n  +  ci££n))A2^  (A44  +  C2(2£>44  +  C2£i4)  +  ck2(£>22  +  ci(2£22 
+  C1H22))  —  (£2  +  Ci(2£4  +  ci£6)V2  —  (£>66  +  ci(2£66  +  ci££66))A2^  ^  —  A55  —  C2(2£>55  +  C2Fqq)  +  ci  ^  —  a2 (£12 
+  2£66  +  ci (££12  +  2Hqq))  +  (£4  +  ciIq)u2  +  (£n  +  ci££n)A2^^/^  —  A*  ^44  +  C2(2£>44  +  C2£44)  +  «2(£>22 
+  Cl  (2F22  +  Cl  ££22))  —  (£2  +  Ci(2£4  +  ciIq))u2  —  (Dqq  +  Ci(2£66  +  CiHqq))X2^  (A55  +  2C2£>55  +  c|£55  +  Cl  ^a2(£i2 
+  2£66  +  Cl  (££12  +  2Hqq))  —  (£4  +  ci£6)oj2  —  (£11  +  ci££n)A2^  +  o?  (£>12  +  Dqq  +  ci(2£i2  +  2£66  +  ci(££i2 
+  H6o)))Xi  (^A44  +  2c2£>44  +  c2£i4  +  ci  [cf2 (F22  +  ^££22)  —  (£4  +  ciIq)uj2  —  (£12  +  2£66  +  Cl  (££12  +  2££66))A2^^ 

with  i  —  1,  •  •  •  ,  12.  When  Eqs.  (67)  are  substituted  into  Eqs.  (66)  a  solution  in  terms 
of  only  twelve  integration  constants  is  obtained.  Thus 

12  12 

Um  (a)  =  ^7 i%  eAiQ  =  ]T nxi  %  eAi“ 

2=1  2=1 

12  12 

Vm  (a)  =  J2Si^i  eAi“  =  2>i  eAl“ 

2=1  2=1 

12  12  12 

wm  (a)  =  J2  Xi  %  eAi Q,  (a)  =  ^  %  eAi  Q,  (a)  =  ^  %  eAi “ 

2=1  2=1  2=1 

(69) 

The  expressions  for  the  generalized  forces  can  then  be  found  in  the  same  way  by  sub¬ 
stituting  Eqs.  (69)  into  Eqs.  (51)  written  in  terms  of  displacements.  Thus,  retaining 
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the  terms  according  to  the  symmetric  cross-ply  composite  shallow  shells  the  following 
expressions  are  derived: 


oo  12 


Nn 


(a,  ft)  =  [An  J2  Xi  7i  Xi  ®i  eAi  a  -  A12  9  E  St  Xi  %  eAi  “ 

I  sin(#  ft)  —  Maa  sin(0  ft) 


m=  1  z=l 

12 


2=1 


+  4r  E  A*  Xi  %  eAi“  +  ttE  a*  ^  ^  eAi< 


2=1 


Npp(a,ft)=Yd 


12 


12 


A66  E  7i  Xi  %  eAi  “  +  A66  E  A*  5*  x*  %  eAi  ‘ 


m=l  L  2=1 

=  A fpp  cos  (9  ft) 


12 


2=1 


12 


cos(#  /?) 


12 


««(“,«=  E  [ffn£  c22  A?  ^  eAi  a  +  2/|  i  E  C2  Af  Xi  ^i  eAi  “  +  i'll  E  c2  X2  %  eXi ‘ 


771=1  2=1 

12 


2=1 


2=1 


12 


12 


2=1 


2=1 


12 


12 


a 


A*  a 


—  2^1 2  E  c2  Aj  <fti  %  eAl  a  —  H12  E  C2  ^2  ^ 2  ^2  eAi  “  -  #12  E  c2  A*  Xi  ^i  eAi 

2=1 
12 

2  2*66  E  C2  A 
2=1 

A55  E  ^  eAl  “  -  A55  E  A*  Ai  eA*  “  —  2  D55  E  C2  ^  e' 


+  4  2266  E  c2  A*  Xi  ^i  eA* 01  ~  2  2^6  E  C2  e  ‘  ”  — 


2=1 


2=1 


„  Ai  01 


2=1 


2=1 


2=1 


12 


12 


12 


2  2A55  E  C2  A,;  Xi  X  eAi  Q  -  F55  E  c2  X  eAi  “  -  ^55  E  c2  A*  Xi  %  eAi  ‘ 


2=1 


x  sin(0  j3)  =  Qa  sin(0  /3) 


12 


2=1 


12 


2=1 


12 


Maa  (a,  /5)  =  E  [Ai  E  A*  eXi  a  +  H^  E  C2  A*  ^  eAl  “  +  E  c2  A?  Xi  @i  e 

772=1  2=1  2=1  2=1 

12  12  12 

+  2  Fn  E c2  A*  ^  eAi“  +  2  Fn  Ec2  X2  Xi  %  eAiQ  -  £>12  E<^  X  eAiQ 


A*  a 


2=1 

12 


2=1 

12 


2=1 


12 


2  F12Ec2  fti  %  eA*Q  -  F12EC2  Xi  ^i  eA‘“  -  #12  E  ci  &  ^  e 


A*  a 


2=1 

12 


2=1 


2=1 


2=1 


H\2  E  c2  Xi  ^i  eAi  “  sin((9  /3)  =  A4aa  sin((9  /3) 
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M, 


PP 


(a,  P)  =  Y  [D™  Y  A*  &  eXi  a  +  D™J2  &  eAi  “  +  Axs  £  &  ^  eAi 

771=1  2=1  2=1  2=1 

12  12  12 

+  D\2  Y^  Aj  fa  3>i  eAl  “  +  2  i^66  ^  C2  3>i  eA*  “  +  2  i<66  ^  C2  fa  3>i  ( 

2=1  2=1  2=1 

12  12  12 

+  2  Fm  Y  c2  A*  x*  %  eA‘  a  +  HmYJ  c2  ^  eAi  “  +  tf66  ^  cl  fa  ^  eA‘ 


Ai  a 


2=1 

12 


2=1 


2=1 


+  2  i?66  ^2  A 


2  X2  ^2  e- 


Ai  a 


2=1 


cos(#  (3)  =  cos(<9  /?) 


(70) 

oo  12  12  12 

Paa  (a,  P)  =  Y  [Hn  Y  ci  Ai  ^  eA*  “  +  Hn  Y  ci  A?  X*  %  eAi  “^n  c2  A,  ^  eAi  “ 

772=1  2=1  2=1  2=1 

12  12  12 

-  Fi2^c2  fa  %  eAi“  -  H12Ycl  <t>i  ®i  eAi“  -  ^i2^ci  Xi  %  eAi" 

2=1  2=1  2=1 

x  sin (0  P)  =  Vaa  sin(0  P) 

Zero  boundary  conditions  are  then  generally  imposed  to  eliminate  the  constants  in  the 
classical  method  in  order  to  establish  the  frequency  equation  for  a  single  plate  element. 
By  contrast,  the  development  of  the  dynamic  stiffness  matrix  entails  imposition  of  gen¬ 
eral  boundary  conditions  in  algebraic  form  which  has  much  wider  implications.  Thus 
in  order  to  develop  the  dynamic  stiffness  matrix  the  following  boundary  conditions  are 
applied: 

Generalized  displacements 

<a  =  0  :  Um  =  Umi  5  Vm  —  Vmi ;  Wm  =  Wmi  5  ^  am  $ai  5  ^  (3  m  ^  (3\  5 

—  W^ml,a 

(71) 

&  —  L  :  Um  —  Um2]  Vm  —  Vm2i  Vi m  —  B/m2 5  ^  am  —  —  ^#27 

W^m,a  =  Wm2a 


Generalized  forces 

Cl  =  0  :  JVaa  =  NqlolY')  ^(3(3  =  ^"(3(3\]  Qa  =  Qa ± 5  *Vlaa  =  •Viaaj.i 

•Viotfi  'M'OifiY  P aa  —  P aa\  5 

(72) 

OL  =  L  \  JVaa  A/*q/q/25  A/^  =  -A/^^2  5  Qa  =  Qq;2  j  JVlaa  A4q,q,2; 

JVla/3  -Vlaf3 2  5  T'aa  P aa2 
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The  imposition  of  Eq.  (71)  leads  to  the  following  relationship: 

8  =  AD 

where 


(73) 


ft  —  {  Umi  Vmi  Wmi  $/9i  ^ml  a  Um 2  Vm2  Wm2  *&a2  *^02  ^m2ta  } 


and 


A  = 


02 

03 

04 

05  ■ 

06  S’-, 

'  08 
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0ii 
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Af 
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Al 

Af 

Af 

Af 
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Af 

Af 
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A 

af 

Af 

Af 

Af 

A 

A 

A 

Aj 

A 

A 

Al 

A 

a? 

Af 

Af 

Af 

-4g 

A 

A 

af 

“4-8 
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Af 

Af 

Af 

A\0 
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Ah 
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Ah 
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Ah 

af0 

a?0 
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ajo 
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Ah 

Ah 

Ah 

Ah 
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Ah 

A  2 

Al  2 

A12 
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af: 2 

ala 

a?2 

a?2 

aj? 

aj-2 

aj? 

(74) 


(75) 


where: 

A\  =  (-l)i+1 

Xi  7i» 

A  =  (-i)i+1  ^  fti, 

A  =  (-!)*+1  Xi 

A\  =  1, 

af  =  (-ir+1  & , 

II 

* 

> 

A  =  (-l)i+1 

Xi  Xi  e(aAi), 

af  =  (-l)i+1  X*  Si  e(“A4, 

af  =  (- l)m  X*  e(“A= 

y. 

O 

II 

> 

aji  =  (-i)i+1  faXi\ 

Af  =  Xi  X,  e(“  A4 

(76) 

with  i  =  1,  •  •  • 

,12.  By  applying  the  same  procedure  for  boundary  conditions  of  the 

generalized  forces  and  thus  exploiting  Eq.  (72),  the  following  relationship  is  obtained: 

F  =  RD  (77) 

where 

F  =  {  A faai  Qai  Aiaai  Aia/g1  ^oai  A/aa2  ^002  Qct2  A4aQ.2  A/la02  Pcm2  } 


£>  =  {  0 1 


04  05  %  07  08  09  01O  0n  012  } J 


(78) 
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and 


n\ 

1Z\ 
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n\ 

n\ 
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n\ 

cj° 

Cl1 
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n\ 

n\ 

n\ 
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H\ o 

Hi  o 

cfo 

Hi  o 

C-Io 

Hi  o 

^?o 

H^ 

H\l 

c 

H\i 

V2 

/cll 

H\  i 

H\ i 

Cfi 

cfi 

n7n 

Cfi 

C?1 

c{? 

V11 

11 

c 

h\2 

V2 

/C12 

h\2 

cf2 

H\2 

C?2 

H"i2 

h\2 

Hi  2 

Cj? 

V11 

7^12 

c 

where 

7^1  =  (_1)Z+1  (fa^l2  O'Xi^i-  A12  Xi  7 i  A*  ~  M2  Xi  “  Ml  Xi  5 
=  (-1)Z+1  ^  -  A66  OLXili-  ^66  Xi  ^  A^  5 


12  “ 
1 

12 

2 

12 

3 

12 

4 
12 

5 

12 

6 

12 

7 
12 

8 

12 

9 

12 

10 
12 
11 
12 
12 


n\  —  —  A55  —  C2  (2  D55  +  C2F55)  —  2  <a2  ci  Fee  —  2  <a2  c2  Hqq  —  Xi  [Ms  F  c2 
(2  F55  +  c2  F55))  A*  -  c?  c2  Xi  Hu  Xi  -  4  a2  c2  Xi#66  A  *  -  a  c\  Fu  fa  A  { 

2  a  ci  FQQcfti  Xi  ci  Ci  H12  (pi  Xi  2  ci  c-^  Fee  (pi  A i  F  ci  Fn  A^  F  <7  i?n  A^ 
+  ci  Xi  TFq  A?; 


(79) 


C4  =  (_1)i+1  («2  C1  Xi  -Cl  2  +  a2  <4  X*  H\2  +  a  C12  ^  +  2aci  Ci2  fa  +  a  c\  H12  fa 
-  Dn  \  ~  2  ci  Fn  Aj  —  Ci  iin  A 4  —  ci  Xi  (-Fn  +  ci  Hu)  \j) ; 

TZ\  —  —a  Dqq  —  2  a  ci  Fqq  —  a  cf  Hqq  —  2  a  ci  Xi  7^66  Xi  —  2  a  cf  Xi  Hqq  Xi  —  (Dqq 
F  ci  (2  Fe6  F  ci  Hqq))  fa  Xp 

K  =  -ci  (-a2  ci  Xi  F12  -  a  7^12  fa  -  a  ci  Hu  fa  F  Fn  A,  F  cx  Fn  A*  F  cx  Xi  Fn  A2); 
Cf  =  (-l)i+1  (  -  A12  a  Xi  e(°  ^  +  An  Xi  e(a  Ai)  7i\  +  A12  Xi  e(a  Ai)  fa 

v  R(3 

+  An  Xi  e(a  Ai)  fa); 

Cf  =  (  — 1)*+1  (^66  OL  Xi  ^  Al^  7 i  +  Xi  $i  c^°  ^  Aj^  ; 

Cf  =  e(“  Ai)  (^55  +  C2  (2  D55  +  c2  F55))  +  2  a2  Cl  e(a  Ai>  C66  +  2  a2  cj  e^a  ^  Hm 

+  X»  e^a  A^  (^55  +  c2  (2  -D55  +  c2  C55))  Aj  +  a2  c2  Xi  e^a  A^  Hi2  \  +  4  a2  c2  %j  <7a  AA 
x  Hqq  Aj  +  a  Cl  e(a  Ai)  Cj2  4>i  \  +  2  a  d  e  L  Xi  F66  fa  Aj  +  a  c2  e(a  A^  H12  fa  Aj 
+  2  a  cj  e(a  Ai)  ii66  fa  Aj  -  Cl  e(“  Ai>  Fn  A2  -  c?  e(a  Ai>  Hn  A2  -  cj  Xi  e(a  Ai)  iin  A?; 
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n{0  =  (-l)i+1  (  -  a2  ci  X*  F\2  -  a2  c2  Xi  e(a  Xi)  H12  -  a  D12  fa  x^  & 

-2  a  cie(a  Xi)  F12  fa  -  a  cj  e(a  Ai)  Hl2  fa  +  Du  e(a  Ai)  +  2  Cl  e(a  Ai)  Fn  fa 
+  cj  e(°  Ai)  Hu  fa  +  ci  x*  e(a  Ai)  (Fn  +  Cl  Hu)  A?) ; 

=  a  Fgg  e(°  Xi) +  2  a  Clfa  Xi )  F66  +  a  c?  e(a  Vtfgg  +  2  a  d  x*  e(a  Ai >  F66  A*  (80) 
+  2  a  Cj  Xi  e*-a  A^  Fgg  A*  +  e^a  x ^  (-Dgg  +  ci  (2  Fgg  +  Ci  Fgg))  q bi  fa ; 

Wj2  =  (-l)i+1  (ci  (-a2  Cl  Xi  e(a  Ai)  H12  -  a  faa  Xi )  F12  fa  -  a  a  e<a  Ai>  Hl2  fa 
+  e(“  Ai)  Fn  Ai  +  ci  e(a  ^FnA*  +  ci  Xi  e<L  A‘>  Fn  A2)) 

with  i  =  1,  •  •  •  ,12.  Now  by  eliminating  the  constants  vectors  D  in  the  Eqs.  (73)  and 
(77)  the  dynamic  stiffness  matrix  which  links  the  forces  and  moments  vector  F  with  the 
generalized  displacements  vector  S  is  derived: 


F  =  K5 ,  K  =  RA1 


fa-l 

fa-l 

faf 

fa f 

fal 

fal 

fal 

faf 

fal 

faf 

faf 

faf 

fa-l 

K2 

fa-l 

fal 

fal 

fal 

fal 

fal 

fa-l 

faf 

faf 

faf 

fa\ 

fa-l 

y 

fal 

fal 

fal 

fa -3 

fal 

fa-l 

faf 

faf 

faf 

fa\ 

fa\ 

fa! 

fa\ 

fal 

fat 

fal 

fa  1 

fa\ 

faf 

faf 

faf 

fa\ 

n 

fal 

fal 

fal 

fal 

fal 

fal 

fa-l 

faf 

faf 

faf 

fa\ 

fa-l 

fa-l 

fal 

fal 

fal 

fal 

fal 

fa-l 

faf 

faf 

faf 

fa\ 

faf 

fal ? 

fal 

fal 

fal 

fal 

fal 

fal 

faf 

faf 

faf 

fa-l 

fa-l 

fal 

fal 

fal 

fal 

fa’s 

fal 

fa-l 

faf 

faf 

faf 

K\ 

fa-l 

fal 

fat 

fal 

fal 

fal 

fa-l 

fa-l 

faf 

faf 

faf 

faf 

faf 

fafo 

fa'll) 

Kfci 

fa% 

fa-10 

fa-io 

fa-io 

fall 

fall 

fall 

fa-h 

fall 

fa-h 

fafl 

fall 

fa'll 

fa\l 

fah 

faf 

fall 

fall 

fa-12 

fa-12 

fa-12 

fa-12 

fa-12 

fa-12 

fa-12 

fa-12 

fa-12 

fall 

fall 

fall 

(81) 


(82) 


The  above  dynamic  stiffness  matrix  will  now  be  used  in  conjunction  with  the  Wittrick- 
Williams  algorithm  [85]  to  analyze  assemblies  of  laminated  composite  cylindrical  and 
spherical  shallow  shells  to  investigate  their  free  vibration  characteristics  based  on  HSDT. 
Explicit  expressions  for  each  shell  element  of  the  DS  matrix  were  obtained  via  symbolic 
computation  using  Mathematica®.  They  are  far  too  extensive  and  voluminous  to  report 
in  this  paper.  The  correctness  of  these  expressions  was  further  checked  by  implement¬ 
ing  them  in  a  MATLAB®  program  and  then  carrying  out  a  wide  range  of  numerical 
simulations. 
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5.4  Assembly  procedure  and  similarities  with  FEM 

Once  the  DS  matrix  of  a  laminated  composite  element  has  been  developed,  it  can  be 
assembled  to  form  the  global  DS  matrix  of  the  final  structure  (see  Fig  6).  Although 


Figure  6:  Direct  assembly  of  dynamic  stiffness  element 

like  the  FEM,  a  mesh  is  required  in  the  DSM,  but  it  should  be  noted  that  the  latter  is 
mesh  independent  in  the  sense  that  additional  elements  are  required  only  when  there  is 
a  change  in  the  geometry  of  the  structure.  A  single  DS  laminate  element  is  enough  to 
compute  any  number  of  its  natural  frequencies  to  any  desired  accuracy,  which,  of  course, 
is  impossible  in  the  FEM.  However,  for  the  type  of  structures  under  consideration  DS 
shell  elements  do  not  have  point  nodes,  but  have  line  nodes  instead.  In  this  particular 
case,  no  change  in  geometry  along  the  longitudinal  direction  is  admitted.  This  is  in 
addition  to  the  assumed  simple  support  boundary  conditions  on  two  opposite  sides, 
inherent  in  DSM  for  shell  elements  at  present.  The  other  two  sides  of  the  shell  can 
have  any  boundary  conditions.  The  application  of  the  boundary  conditions  of  the  global 
dynamic  stiffness  matrix  involves  the  use  of  the  so-called  penalty  method.  This  consists 
of  adding  a  large  stiffness  to  the  appropriate  leading  diagonal  term  which  corresponds 
to  the  degree  of  freedom  of  the  node  that  needs  to  be  suppressed.  It  is  thus  possible  to 
apply  free  (F),  simple  support  (S)  and  clamped  (C)  boundary  conditions  on  the  structure 
by  penalizing  the  appropriate  degrees  of  freedom.  Clearly  for  simple  support  boundary 
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condition,  the  generalized  displacement  amplitudes  Vi,  Wi  and  are  assigned  zero 
values.  On  the  other  hand,  for  clamped  boundary  condition  Ui,  Vi,  Wi,  <f>Xi,  <hyi  and 
WXi  on  the  boundary  are  assigned  zero  values.  Of  course  for  the  free-edge  boundary 
condition  stress  resultants  are  assigned  zero  values  and  then  no  penalty  will  be  applied 
at  the  generalized  displacement  amplitudes.  Because  of  the  similarities  between  DSM 
and  FEM,  DS  elements  can  be  implemented  in  FEM  codes  to  enhance  the  accuracy  of 
results  in  FEM  very  considerably. 

6  Results  and  discussion 

Free  vibration  characteristics  of  laminated  composite  cylindrical  and  spherical  shallow 
shells  with  symmetric  cross-ply  lamination  schemes  are  investigated.  Results  have  been 
obtained  in  exact  sense  by  using  the  Levy-type  closed  form  solution  within  the  framework 
of  the  DSM.  The  material  used  in  all  the  analysis  carried  out  is 

=  25;  G12  =  Gis  =  0.5  E2 ;  G23  =  0.2  E2]  v12  =  0.25; 

E/2 

if  not  differently  stated.  A  first  validation  is  undertaken  in  Table  4,  the  fundamental 
circular  frequency  parameter  of  three-layer  symmetric  cross-ply  and  moderately  thick 
spherical  and  cylindrical  shallow  shells  has  been  computed.  The  results  obtained  by 
using  the  present  DSM  formulation  have  been  compared  with  those  proposed  by  Khdeir 
and  Reddy  [87]  by  using  the  classical  Levy-type  solution.  As  can  be  seen  in  Table  4  the 

2  / 

Table  4:  Fundamental  circular  frequency  parameter  Co  =  uo  w  of  square  cylindrical 

and  spherical  shells  with  staking  sequence  [0°/90°/0°]  length-to-thickness  ratio  a/h  =  10 
and  radius-to-length  ratio  Rp/a  —  20. 


ssss 

scss 

SFSF 

scsc 

SFSS 

SFSC 

Spherical  shell 

Khdeir  and  Reddy  [87] 
Present  HSDT  DSM 

11.807 

11.807 

13.481 

13.480 

3.797 

3.796 

16.100 

16.100 

4.328 

4.328 

6.088 

6.088 

Cylindrical  shell 

Khdeir  and  Reddy  [87] 
Present  HSDT  DSM 

11.793 

11.793 

13.825 

13.825 

3.789 

3.789 

15.999 

15.999 

4.322 

4.322 

6.089 

6.089 

results  perfectly  match  each  other  for  all  the  considered  boundary  conditions.  In  Tables 
5  and  6  further  assessments  have  been  carried  out.  Most  notably,  the  fundamental  circu¬ 
lar  frequency  parameter  has  been  computed  for  three  and  four-layer  symmetric  cross-ply 
and  moderately  thick  spherical  and  cylindrical  shallow  shells,  respectively.  Results  are 
compare  with  several  meshless  formulations  proposed  by  Ferreira  el.  all  [68,88,89],  in 
particular  a  radial  basis  formulation  (RBF)  based  on  a  HSDT,  a  wavelet  collocation 
(WLC)  formulation  based  on  a  FSDT  and  a  radial  basis  collocation  (RBFC)  solution 
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2  / 

Table  5:  Dimensionless  fundamental  circular  frequency  parameter  Co  —  of 

square  spherical  and  cylindrical  shells  with  staking  sequence  [0°/90°/0°],  length-to- 
thickness  ratio  a/h  =  10  and  varying  the  radius-to-length  ratio  Rp/a  and  R/a. 


Shell  Configuration 

Theory 

R/a 

5 

10 

20 

50 

100 

Plate 

Spherical 

HSDT  RBF  [68] 

12.063 

11.861 

11.810 

11.796 

11.794 

11.793 

FSDT  WLC  [88] 

12.417 

12.227 

12.179 

12.165 

12.164 

12.163 

SSDT  RBFC  [89] 

12.125 

11.966 

11.925 

11.914 

11.912 

11.911 

HSDT  DSM 

12.054 

11.857 

11.807 

11.793 

11.790 

11.790 

R e/ 

'a 

5 

10 

20 

50 

100 

Plate 

Cylindrical 

HSDT  RBF  [68] 

11.851 

11.808 

11.797 

11.794 

11.793 

11.793 

FSDT  WLC  [88] 

12.214 

12.176 

12.166 

12.163 

12.163 

12.163 

SSDT  RBFC  [89] 

11.923 

11.915 

11.913 

11.912 

11.912 

11.912 

HSDT  DSM 

11.846 

11.804 

11.793 

11.790 

11.790 

11.790 

2  / 

Table  6:  Dimensionless  fundamental  circular  frequency  parameter  Co  =  ojy  of 

square  spherical  and  cylindrical  shells  with  staking  sequence  [0°/90o/90o/0°],  length-to- 
thickness  ratio  a/h  —  10  and  varying  the  radius-to-length  ratio  Rp/a  and  R/a. 


Shell  Configuration 

Spherical 

Theory 

HSDT  RBF  [68] 
FSDT  WLC  [88] 
SSDT  RBFC  [89] 
HSDT  DSM 

R/a 

5 

12.051 

12.493 

12.099 

12.043 

10 

11.848 

12.299 

11.938 

11.844 

20 

11.796 

12.250 

11.896 

11.793 

50 

11.782 

12.236 

11.885 

11.779 

100 

11.780 

12.234 

11.883 

11.777 

Plate 

11.779 

12.233 

11.883 

11.776 

Rp/a 

5 

10 

20 

50 

100 

Plate 

Cylindrical 

HSDT  RBF  [68] 

11.838 

11.794 

11.783 

11.780 

11.779 

11.779 

FSDT  WLC  [88] 

12.279 

12.240 

12.230 

12.228 

12.227 

12.227 

SSDT  RBFC  [89] 

11.901 

11.887 

11.884 

11.883 

11.883 

11.883 

HSDT  DSM 

11.832 

11.790 

11.780 

11.777 

11.777 

11.776 
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procedure  based  on  a  sinusoidal  shear  deformation  theory  (SSDT).  As  can  be  seen  in 
Tables  5  and  6  the  present  formulation  leads  to  results  which  are  in  an  excellent  agree¬ 
ment  with  those  obtained  by  using  the  aforementioned  theories,  for  different  values  of 
the  radius-to-length  ratio  ( R/a ).  From  Tables  7  to  10  the  investigation  is  focused  on 
symmetric  cross-ply  cylindrical  shells,  two  different  stacking  sequences  are  investigate, 
and  both  moderately  thick  and  thin  shallow  shells  are  taken  into  account.  The  first  five 
circular  frequency  parameters  are  computed  for  several  boundary  conditions.  From 
Fig.  7  to  Fig.  9  the  first  six  mode  shapes  of  SCSC,  SFSC  and  SFSF  square  cylindrical 
shallow  shells  are  depicted.  Most  notably,  the  mode  shapes  are  relative  to  a  symmetric 
cross-ply  [0°/90°/0°]  stacking  sequence,  a/h  =  100  and  R/a  —  30.  An  assessment  of 
the  present  formulation  is  carried  out  in  Table  6  for  isotropic  spherical  panels.  The 


2  / 

Table  7:  First  five  circular  frequency  parameters  Co  —  u  . N/r,  °f  square  cylindrical 

shells  with  staking  sequence  [0°/90°/0°]  length-to-thickness  ratio  a/h  —  10  and  varying 
the  radius-to-length  ratio  Rp/a. 


Mode 

SSSS 

SFSS 

SFSF 

Ftp/ a 

5 

50 

100 

5 

50 

100 

5 

50 

100 

1 

11.846 

11.790 

11.790 

4.376 

4.326 

4.326 

3.783 

3.790 

3.790 

2 

18.489 

18.487 

18.487 

14.398 

14.401 

14.401 

5.726 

5.576 

5.574 

3 

29.399 

29.333 

29.332 

17.307 

17.219 

17.218 

13.898 

13.905 

13.905 

4 

30.856 

30.860 

30.860 

22.971 

22.959 

22.959 

15.741 

15.734 

15.734 

5 

32.961 

32.949 

32.948 

28.330 

28.336 

28.336 

23.019 

23.020 

23.020 

Mode 

SFSC 

SCSS 

SCSC 

Ftp /a 

5 

50 

100 

5 

50 

100 

5 

50 

100 

1 

6.129 

6.090 

6.089 

13.865 

13.823 

13.823 

16.026 

15.997 

15.997 

2 

14.978 

14.981 

14.981 

19.722 

19.720 

19.720 

21.131 

21.130 

21.130 

3 

18.869 

18.793 

18.793 

30.980 

30.920 

30.920 

32.306 

32.310 

32.310 

4 

24.118 

24.107 

24.107 

31.524 

31.528 

31.528 

32.467 

32.414 

32.413 

5 

28.594 

28.600 

28.600 

34.318 

34.308 

34.308 

35.629 

35.618 

35.618 

first  six  natural  frequencies  are  computed  and  compared  with  the  3D  elasticity  solution 
and  many  other  theories,  which  are  refined  in  the  displacement  field  or  in  the  curvature 
description.  As  can  be  seen  from  the  Table  the  proposed  DSM  formulation  provide  the 
best  accuracy  for  the  fundamental  natural  frequency  with  respect  to  the  3D  elasticity 
solution,  but  for  higher  frequencies  there  is  an  increase  in  the  error  percentage.  This 
slight  loss  of  accuracy  is  due  to  the  fact  that  approximated  curvature  descriptions  have 
been  employed.  In  Table  12  the  first  three  circular  frequency  parameters  of  three-layer 
symmetric  cross-ply  shallow  spherical  shells  are  calculated.  The  boundary  conditions 
SSSS,  SCSS  and  SCSC  are  taken  into  account.  The  investigation  is  carried  out  for  dif¬ 
ferent  values  of  the  radius-to-length  ratio  and  for  moderately  thick  and  thin  spherical 
panels.  In  Table  13  the  effect  of  the  orthotropic  ratio  on  the  dimensionless  fundamen¬ 
tal  frequency  parameter  is  studied.  The  boundary  conditions  SSSS,  SCSS  and  SCSC 


41 


Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


FA8655-10-1-3084  Report  6  Dynamic  Stiffness  Modelling  of  Plate  and  Shell  Assemblies 


(a)  Mode  1,  w  =  31.916 


(b)  Mode  2,  Cj  =  36.277 


(c)  Mode  3,  w  =  49.256 


(d)  Mode  4,  w  =  72.406 


(e)  Mode  5,  Cj  =  84.982  (f)  Mode  6,  w  =  87.738 

Figure  7:  First  six  mode  shapes  of  a  symmetric  cross-ply  cylindrical  shell  with  SCSC  boundary 
condition. 
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(a)  Mode  1,  Cj  =  7.031 


(c)  Mode  3,  Cj  —  32.260 


(b)  Mode  2 ,  u  =  17.257 


(d)  Mode  4,  Cj  —  36.340 


(e)  Mode  5,  c5  =  37.804  (f)  Mode  6,  c j  =  51.898 

Figure  8:  First  six  mode  shapes  of  a  symmetric  cross-ply  cylindrical  shell  with  SFSC  boundary 
condition. 
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(a)  Mode  1,  Cj  =  3.919 


(c)  Mode  3,  w  =  15.649 


(b)  Mode  2 ,  u>  =  6.334 


(d)  Mode  4,  Ci  =  18.007 


(e)  Mode  5,  Cj  =  33.435  (f)  Mode  6,  w  =  35.162 

Figure  9:  First  six  mode  shapes  of  a  symmetric  cross-ply  cylindrical  shell  with  SFSF  boundary 
condition. 
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Table  8:  First  five  circular  frequency  parameters  Co  =  uo  ^  of  square  cylindrical 

shells  with  staking  sequence  [0°/90°/0°]  length-to-thickness  ratio  a/h  =  100  and  varying 
the  radius-to-length  ratio  Rp/a. 


Mode 

ssss 

SFSS 

SFSF 

Rp/a 

5 

50 

100 

5 

50 

100 

5 

50 

100 

1 

20.329 

15.238 

15.193 

8.271 

4.556 

4.511 

3.913 

3.921 

3.913 

2 

23.745 

22.818 

22.811 

16.599 

16.270 

16.262 

10.742 

6.061 

5.950 

3 

40.355 

56.135 

40.153 

31.030 

23.301 

23.231 

15.657 

15.649 

15.649 

4 

61.528 

60.093 

56.089 

32.099 

30.232 

30.218 

19.269 

17.984 

17.976 

5 

61.560 

66.369 

60.084 

35.833 

35.773 

35.773 

35.160 

33.247 

33.173 

Mode 

SFSC 

scss 

scsc 

Rp/a 

5 

50 

100 

5 

50 

100 

5 

50 

100 

1 

10.746 

6.933 

6.888 

26.683 

22.586 

22.549 

34.921 

31.856 

31.833 

2 

17.765 

17.249 

17.249 

29.277 

28.410 

28.402 

36.984 

36.266 

36.259 

3 

36.438 

32.124 

32.071 

43.951 

43.726 

43.726 

49.460 

49.251 

49.251 

4 

38.696 

36.341 

36.341 

68.784 

68.716 

68.716 

72.476 

72.409 

72.409 

5 

39.498 

37.779 

37.762 

74.430 

69.904 

69.866 

88.755 

84.909 

84.880 
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Table  9:  First  five  circular  frequency  parameters  Co  =  uo  ^  of  square  cylindrical 

shells  with  staking  sequence  [0°/90o/90o/0°]  length-to-thickness  ratio  a/h  =  10  and 
varying  the  radius-to-length  ratio  Rp/a. 


Mode 

SSSS 

SFSS 

SFSF 

Rp/a 

5 

50 

100 

5 

50 

100 

5 

50 

100 

1 

11.832 

11.777 

11.777 

5.782 

5.750 

5.749 

5.367 

5.377 

5.378 

2 

21.836 

21.837 

21.837 

16.726 

16.631 

16.631 

6.826 

6.700 

6.699 

3 

27.461 

27.381 

27.381 

18.984 

18.991 

18.991 

18.659 

18.668 

18.668 

4 

33.244 

33.232 

33.232 

25.249 

25.240 

25.240 

19.895 

19.893 

19.893 

5 

37.445 

37.452 

37.452 

33.637 

33.538 

33.538 

22.875 

22.870 

22.869 

Mode 

SFSC 

SCSS 

SCSC 

Rp/a 

5 

50 

100 

5 

50 

100 

5 

50 

100 

1 

7.049 

7.018 

7.018 

13.491 

13.446 

13.446 

15.304 

15.271 

15.270 

2 

17.991 

17.906 

17.905 

22.691 

22.693 

22.693 

23.694 

23.695 

23.695 

3 

19.371 

19.377 

19.377 

28.772 

28.699 

28.698 

30.001 

29.934 

29.933 

4 

26.068 

26.060 

26.060 

34.292 

34.281 

34.281 

35.304 

35.294 

35.294 

5 

35.144 

35.053 

35.053 

37.899 

37.905 

37.905 

38.435 

38.442 

38.442 
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Table  10:  First  five  circular  frequency  parameters  Co  =  uo  ^  w  of  square  cylindrical 

shells  with  staking  sequence  [0°/90o/90o/0°]  length-to-thickness  ratio  a/h  —  100  and 
varying  the  radius-to- length  ratio  R^/a. 


Mode 

ssss 

SFSS 

SFSF 

Rp/a 

5 

50 

100 

5 

50 

100 

5 

50 

100 

1 

20.359 

15.231 

15.186 

9.198 

6.156 

6.126 

5.692 

5.700 

5.700 

2 

28.604 

27.864 

27.857 

23.409 

22.766 

22.691 

15.571 

7.353 

7.253 

3 

54.573 

53.795 

53.748 

30.923 

23.184 

23.184 

22.758 

22.758 

22.758 

4 

59.700 

54.436 

54.439 

35.519 

33.831 

33.820 

25.242 

24.420 

24.409 

5 

61.548 

60.065 

60.056 

51.501 

51.471 

51.471 

41.995 

32.166 

32.082 

Mode 

SFSC 

scss 

SCSC 

Rp/a 

5 

50 

100 

5 

50 

100 

5 

50 

100 

1 

11.411 

7.942 

7.898 

26.272 

22.026 

21.996 

33.964 

30.720 

30.697 

2 

24.209 

23.835 

23.835 

32.962 

32.215 

32.207 

39.376 

38.726 

38.718 

3 

38.139 

31.045 

30.985 

57.062 

56.905 

56.905 

61.024 

60.867 

60.867 

4 

41.815 

40.169 

40.157 

71.834 

66.814 

66.777 

85.287 

80.995 

80.958 

5 

51.904 

51.837 

51.837 

73.372 

72.054 

72.045 

86.582 

85.426 

85.419 

are  accounted  for,  a  symmetric  cross-ply  [0°/90o/90o/0°]  square  spherical  shallow  shell 
(R/a  =  50)  is  analyzed.  As  expected  the  fundamental  frequency  parameter  increases 
when  increasing  both  the  orthotropic  ratio  and  the  length-to-thickness  ratio,  for  all  the 
considered  boundary  conditions.  Finally,  in  Fig.  10  the  first  six  mode  shapes  of  SCSC 
square  spherical  shallow  shells  are  represented.  The  geometrical  characteristics  and  the 
lamination  scheme  are  equal  to  the  ones  used  for  the  cylindrical  shallow  shells  shown  in 
Figs.  7,  8  and  9. 

7  Conclusions 

A  detailed  literature  reviews  on  the  dynamic  stiffness  method  for  both  buckling  of  lam¬ 
inated  composite  plates  and  free  vibration  of  doubly  curved  shallow  shells  have  been 
presented.  In  particular,  in  the  first  part  of  the  report  the  stability  equations  have  been 
obtained  using  the  principle  of  minimum  potential  energy  and  the  dynamic  stiffness 
matrix  has  been  derived  for  laminated  composite  plate  elements  based  on  the  HSDT. 
The  element  stiffnesses  have  been  implemented  in  a  computer  program  and  results  for 
composite  plate  assemblies  have  been  obtained  and  validated. 

In  the  second  part,  an  exact  free  vibration  analysis  of  laminated  composite  shallow 
shells  has  been  carried  out  by  combining  for  the  first  time  the  dynamic  stiffness  method 
and  a  higher  order  shear  deformation  theory.  The  effect  of  several  parameters  such  as 
length-to-thickness  ratio  and  radius-to-length  ratio,  orthotropic  ratio,  stacking  sequence 
and  number  of  layers  on  the  dimensionless  circular  frequency  parameters  has  been  in¬ 
vestigated  in  details.  Results  have  been  compared  with  those  available  in  the  literature 
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Table  11:  First  six  natural  frequencies  Umn  rad/ s,  of  square  isotropic  spherical  shells. 

Natural  frequencies 


Theory 

uji  (1,1) 

w2  (2,1) 

u,3  (1,2) 

0,4  (2,  2) 

u,5  (3,1) 

o,6  (1,3) 

3D  [90] 

0.52543 

0.58420 

0.58487 

0.67676 

0.75219 

0.75220 

47 

CST  [91] 

FSTD  [92] 

HSTD  [92] 

FSTD  (DT*)  [93] 
FSTD  (ST*  )  [93] 

0.53263 

0.50211 

0.50223 

0.52864 

0.52830 

(-1.370)t 

(+4.438) 

(+4.415) 

(-0.611) 

(-0.546) 

0.59041 

0.56247 

0.56276 

0.58954 

0.58853 

(-1.063) 

(+3.720) 

(+3.670) 

(-0.914) 

(-0.741) 

0.59080 

0.56248 

0.56277 

0.58954 

0.58853 

(-1.014) 

(+3.828) 

(+3.779) 

(-0.798) 

(-0.626) 

0.68486 

0.65706 

0.65788 

0.68370 

0.68232 

(-1.197) 

(+2.911) 

(+2.790) 

(-1.025) 

(-0.822) 

0.76020 

0.73915 

0.73966 

0.75974 

0.75818 

(-1.065) 

(  +  1.734) 

(  +  1.666) 
(-1.004) 
(-0.796) 

0.76260 

0.74035 

0.74081 

0.75974 

0.75818 

(-1.383) 

(  +  1.575) 

(  +  1.514) 
(-1.002) 
(-0.795) 

HSDT  DSM 

0.52795 

(-0.480) 

0.58899 

(-0.820) 

0.58982 

(-0.846) 

0.68562 

(-1.309) 

0.75989 

(-1.023) 

0.76032 

(-1.080) 

Material  and  geometric  properties:  a  =  b  =  1.0118,  h  =  0.0191,  R  —  1.91,  E  —  1,  p  —  1,  v  —  0.3 
f  Error%  =  U  X  100;  |  Donnell  approximation;  *  Sanders  approximation 
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Table  12:  First  three  circular  frequency  parameters  Co  =  oo  ^  of  square  spherical 

shells  with  staking  sequence  [0°/90°/0°]  and  varying  the  radius-to-length  R/a  and  the 
length-to-thickness  a/h  ratios. 


Mode 

a/h 

ssss 

scss 

scsc 

10 

Rp/a 

20 

50 

100 

20 

50 

100 

20 

50 

100 

1 

8.611 

8.601 

8.600 

9.999 

9.990 

9.988 

11.494 

11.485 

11.484 

2 

14.572 

14.563 

14.562 

15.345 

15.335 

15.335 

16.237 

16.227 

16.226 

3 

20.876 

20.868 

20.868 

21.989 

21.982 

21.981 

23.038 

23.031 

23.030 

Mode 

SSSS 

SCSS 

SCSC 

50 

Rp/a 

20 

50 

100 

20 

50 

100 

20 

50 

100 

1 

11.171 

10.951 

10.923 

15.972 

15.766 

15.735 

21.830 

21.649 

21.624 

2 

18.277 

18.078 

18.050 

21.628 

21.419 

21.390 

26.255 

26.049 

26.017 

3 

33.506 

33.336 

33.311 

35.499 

35.329 

35.304 

38.481 

38.297 

38.272 

Mode 

SSSS 

SCSS 

SCSC 

100 

R/3  /  a 

20 

50 

100 

20 

50 

100 

20 

50 

100 

1 

10.279 

11.187 

11.067 

14.130 

16.314 

16.201 

23.562 

22.824 

22.718 

2 

16.662 

18.357 

18.243 

22.832 

21.981 

21.860 

28.044 

27.207 

27.087 

3 

34.561 

33.895 

33.795 

36.767 

36.086 

35.987 

40.249 

39.525 

39.419 

2  / 

Table  13:  Fundamental  circular  frequency  parameter  Co  —  oo  ^  \/e^i  °f  square  spherical 

shells  with  staking  sequence  [0°/90o/90o/0°],  R/a  =  50  and  varying  orthotropic  E1/E2 
and  the  length-to-thickness  a/h  ratios. 


BCs 

a/h 

Ei/E2 

5 

10 

20 

30 

40 

SSSS 

5 

6.4534 

7.212 

8.038 

8.530 

8.879 

10 

7.697 

9.201 

11.110 

12.336 

13.223 

50 

8.410 

10.511 

13.738 

16.288 

18.445 

100 

8.702 

10.803 

14.050 

16.660 

18.904 

sssc 

5 

7.134 

7.812 

8.563 

9.047 

9.414 

10 

9.427 

11.102 

12.877 

13.915 

14.653 

50 

11.207 

14.585 

19.467 

23.140 

26.136 

100 

11.495 

14.957 

20.098 

24.101 

27.467 

scsc 

5 

7.881 

8.525 

9.254 

9.752 

10.148 

10 

11.360 

13.112 

14.764 

15.686 

16.341 

50 

14.903 

19.817 

26.590 

31.460 

35.296 

100 

15.259 

20.432 

27.881 

33.563 

38.282 
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(a)  Mode  1,  u  =  32.361 


(b)  Mode  2,  u  =  36.905 


(c)  Mode  3,  =  49.861  (d)  Mode  4,  =  72.907 


(e)  Mode  5,  Co  =  85.063  (f)  Mode  6,  Cj  —  87.862 

Figure  10:  First  six  mode  shapes  of  a  symmetric  cross-ply  spherical  shell  with  SCSC  boundary 
condition. 
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including  3D  elasticity  solution.  Different  boundary  conditions  have  been  considered, 
according  to  the  Levy-type  closed-form  solution.  Representative  mode  shapes  have  been 
presented  and  discusses.  The  investigation  has  shown  that  the  DSM  allows  computing  all 
the  natural  frequencies  of  thin  and  thick  laminated  composite  cylindrical  and  spherical 
shallow  shells  with  high  accuracy  and  low  computational  cost.  The  implementation  of 
the  HSDT  within  the  framework  of  the  DSM  permits  considerable  improvements  in  re¬ 
sults  accuracy  when  compared  to  other  analytical,  FEM  or  meshless  formulations  based 
on  FSDT. 

All  the  Work  packages  have  been  successfully  completed  (see  Gantt  chart  in  Fig.  11  and 
the  proposal  [94]).  The  project  is  completed  as  scheduled. 
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Appendix  A:  Laminate  geometric  and  constitutive  equa¬ 
tions 


The  geometric  relation  for  a  lamina  in  the  local  or  lamina  reference  system  can  be 


£xx 

T)x  0  0 

£yy 

0  T>y  0 

u 

Ixy 

= 

T>y  T>x  0 

V 

lyz 

0  T)Z  T>y 

w 

Ixz  - 

T>z  0  T>x 

and  in  terms  of  the  functional  degrees  of  freedom: 


written 


(A.l) 


£xx 

£yy 

Ixy 

lyz 


T>x  0  fci  zA  T)xx 

0  Tci  Z^J  ^yy 

r^y  Vx  fci  z  r>Xy  -\-  c\  z  r)yX^ 

0  0  (l  +  3  Cl  z2\  T>y 

0  0  (l  +  3  ci  z2)  T>x 


(*  +  ci  z3)  T>x 
0 

(z  +  ci  z3)  T>y 
0 

(l  +  3  ci  z2) 


0 

z  +  c  1  z3)  T>y 
z  +  ci  z3J  T>x 
1  +  3  ci  z2) 

0 


(A- 2) 


where  Vx  and  Vy  are  the  derivatives  in  x  and  y  respectively  and  ci  =  —yjy-  The 
constitutive  equations  in  the  lamina  reference  system  can  be  written,  in  terms  of  reduced 
stiffness  coefficients,  as: 


cri 

On 

012 

0 

0 

0 

C2 

Cl  2 

C22 

0 

0 

0 

T1 2 

= 

0 

0 

066 

0 

0 

T~23 

0 

0 

0 

044 

0 

-  T13  - 

0 

0 

0 

0 

055 

(A.3) 


where  the  Cij  are  expressed  in  terms  of  stiffness  coefficients  CV/,  as: 


r2 

O13 


E1 


Cn  =  Cn  -  ^  , 

O33  1  —  ^12^21 


fy  fy  C\3  C23  VI2E2  ~  C23  E2 

^12  ^12  T  5  ^22  ^22  Tj  T 

O33  1  —  ^12^21  O33  1  —  ^12^21 


C44  =  C44  =  G  23,  C55  —  C55  =  tr  13  C66  —  C*66  =  G 12 

(A. 4) 

where  £1  is  the  elastic  modulus  in  the  fibre  direction,  i?2  the  elastic  modulus  in  per¬ 
pendicular  to  the  fibre,  v\2  and  z/21  =  V12E2/E1  the  Poisson’s  ratios,  G 12  =  G13  and 
G23  the  shear  modulus  of  each  single  orthotropic  lamina.  If  the  lamina  is  placed  at  an 
angle  6  in  the  laminate  or  global  reference  system,  the  equation  need  to  be  transformed 
as  follows: 


Cn  =CnC4  +  2(Ci2  +  2  Cm)S2C2  +  C22S 4 
Ci2  =(Cn  +  C22  —  4C66)*S2C2  +  Ci2(<S4  +  C4) 

Ci6  =(Cn  -  Cn  -  2 C66)SC3  +  (C12  -  C22  +  2C66)«S3 
C22  =CnS4  +  2(Ci2  +  2C66)52C2  +  C22C4 

C26  =(Cn  -  C12  -  2C66)53C  +  (C12  -  C22  +  2C66)SC3  (A. 5) 

C66  =(Cn  +  C22  -  2Ci2  -  2C66)<S2C2  +  C66(<S4  +  C4) 

C44  =C44C2  +  C55S2 
C55  =C44<S2  +  C55C2 
C45  =(C55  -  CU)CS 
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where  C  =  cos  ( 9 )  and  S  —  sin  (6).  This  leads  to  the  constitutive  equation  for  the  k- th 
lamina  in  the  laminate  or  global  reference  system: 

o XX  r  c u  c  12  c is  o  o 

CJyy  C 12  C  22  C*26  0  0 

TXy  —  C\Q  C  26  £*66  ^  0 

Tyz  0  0  0  C44  C45 

T’sz  _  0  0  0  C*45  C^55 

that  in  compact  form  can  be  written  for  each  k- th  lamina  as: 

£ Tk  =  Ck£k  (A.7) 
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Appendix  B:  Polynomials  Coefficients 

The  polynomial  coefficients  are  following  defined: 

0-1  =  c\  {F2X  -  £>n  iJn)(£>66  +  Cl  (2  F66  +  Cl  H66 )) 

a2  =  -cl  (a44  +  a2£>22)  F^  +  2a2c^i2FnFi2  -  c?  c\dxxf12  +  2a2c?FnF122  -  2a2c?F21F22  +  4a2c2£>i2 

FhF66  -  4a2c2£>nFi2F66  -  4a2c2£>nFg6  -  8a2cfFnFg6  +  A44c2£>nFfn  -  a2c2£>22FJn  +  a2c2£>n £>22Ffn 
-  2cc2  c\d12DqqHxx  -  2a2c?£>i2Fi2FJn  -  4a2c?£>66Fi2FJii  -  a^F^FFn  +  2a2c?DiiF22Hn  -  4a2c4Fi2F66 
Hu  -  4a2  c4  Fqq  f/ 1 1  +  2cc2  c\dxxDqqH12  +  2a2c?£>i2Fn  FJi2  +  4a2c?%Fn  Ffi2  -  2a2cf  D11F12H12  +  2a2c4Fn 
Fi2FJi2  +  4cx2  cl  F\\FqqH\2  —  a2  c4£>nFJ22  —  a2c^F21F22  +  a2c4£>n  FJh  FF22  +  4a2  c2£>n  £>66^66  +  4a2cf  Di2Fn 
#66  +  8a2cf  £>66FhFF66  —  4a2c^  £>n  Fi2FF66  +  2a2c4£>i2FJnFF66  +  4a2c4£>66#ll#66  —  2a2c4£>nFJi2FJ66  +  A55£>n 
(#66  +  ci(2F66  +  ci Hqq))  +  c2  ^£>h£>66F55  +  2ci£>iiF55F66  +  c2  (  —  F2-LF44  +  £>nF44FJn  +  £>iiF55FF66^  +  2c2 
^-c2£>44F21  +  £>n  (c2£>44FFh  +  £>55  ^£>66  +  2ciF66  +  c2#66)))  +  (#11  +  C]  (2Fn  +  Ci  Ffn  ))  (£>66  +  ci  (2F66  +  ciHqq))XNx 


-  z,uc  c2  J^l 


a-3  =  2a  c2 £>12 £>55  _  _ _ _  ^  _  _  _  _  _  _ _  _  _  _ _ 

—  8a2cic2£>44£>66Fn  +  4a2cic2£>n£>44£i2  +  4a2cic2£>i2£>55£i2  +  8a2cic2£>55£>66Fi2  +  2a4c2F22^n^i2  , 
c2c2£>44£ii£i2  -  2a4c2£>i2£lf2  +  2a2c2c2£>55F22  -  4 ol4  c\f^2  -  4o2c1c2D11Dq5F22  -  2a4c2£>i2FnF22  -  2 
cdc2  £>66  Fi  i  F22  +  2a4c2£>nFi2F22  +  4cx4c\f11F12F22  —  2c2£>n£>55F44  —  a2c2£>n  £>66#44  —  2a2cic2£>i2Fn 
F44  —  4a2cic2£>66FnF44  +  2a2cic2£>nFi2F44  +  2a2c2c2FnFi2F44  +  cx2  c^D^Fqq  —  a2c2£>n£>22£55  —  2c2 
#11#44#55  +  2a2c2£>i2£>66F55  +  2a2cic2£>i2Fi2F55  +  4a2cic2£>66#l2F55  +  a2c2c2F22F5  5  —  2a2cic2£>n 
F22F55  —  c2£>iiF44F55  +  4a2cic2FnF44F66  +  4a4c2£>22Fn  F66  +  8a2c2c2£>44FnF66  —  8a4c2£>i2Fi2F66  +  8a2c2 
c2#55#l2#66  —  16a4cfF22F66  |4a4c2DnF22F66  +  12a4cf FnF22F66  +  2a2cic2£>nF44F66  +4 a2c2c2FnF44 

F66  +  4a2c2c2Fi2F55F66  -  8a4c2£>i2F26  +  8a2c2c2F55F26  -  16a4c?Fi2F26  +  4a2c2c2F55F26  -  4a2c2 
c2#i2#44#ll  —  cx4c\d22DqqHh  —  8a2c2c2F44Fe6Fii  +  2cx4c\d22F\2H\\  —  2cx4c\D\2F22H\\  —  4cx4c\DqqF22 
#11  +  2a4  c4  F\2F22H\\  —  2a2c2c2Fi2F44FFn  —  4a2c2c2£>66F44FJn  +  2o4c\d22FqqH11  +  4a4c4F22F66FFn  +  2 
a4c2F22FFi2  —  2a4c2£>n£>22FJi2  +  4cx4c\d12DqqHi2  —  2o.4c\d22F11Hi2  +  4a4c^  £>i2Fi2FJi2  +  8a4c^F66£l2 
H12  -  2a4c4F22Fi2  -  2a4cf  DxlF22Hr%  -  2a4c4Fn F22 ffi2  -  8a4c4Fi2F66Fi2  -  8a4c4F26Fi2  +  2a4c4£>i2 
H212  +  4a4c4D66F22  —  2a2c2c2£>n£>55FF22  —  a4c2£>n£>66#22  —  2a4c^Fi2FnH22  —  4a4Ci£>66#ll#22  +  2  a4cf 
DhF12H22  +  2a4  c4  F\\F\2H22  —  cx2c\c^DhFqqH22  +  2a4cf  FnF66F22  +  4a4c4FnFe6#22  —  2a4c4£>i2FJnFJ22 

—  4a4c4£>66#ll#22  +  4a4c2£>22  FJ66  —  4a4c2DnD22F66  —  2  a2c2c2£>n  £>44Ff66  —  4a2c2c2£>i2£>55FF66  +  8a4c2£>i2 
#66#66  —  8a2  c1c2D55D66Hqq  —  4(x4c\d22Fi\Hqq  +  8a4cf  Fi2Fi2FF66  +  16a4cf  DqqF12Hqq  —  4  ol4c\D\iF22Hqq 

—  2a4c4FnF22FF66  —  0'2c2c2FnF44FF66  —  2a2c2c2£>i2F55FF66  —  4a2c2c2F66F55FF66  —  a4  c4F22FFnFF66  +  4a4c4 
D\2H\2Hqq  +  8a4c4£>66-FFi2F66  —  ck4c4FiiFF22FF66  +  cx2  d\2\N  xq  —  ex2  D\\D22\Nxq  —  2c2Di\D^\Nxq  +  2cx2  D\2Dqq 
^NxO  ~  2c2D55DqqXNx0  —  2cx2  c1D22F1i  XNx0  —  4ci  c2F44Fn  A  +  4a2ci£>i2Fi2AiV£Co  +  4a2ciF66Fi2A 

Nxo  +  4a2c2F-j22AiV£Co  —  2a2ciDnF22ANa;o  —  4a2 c2FhF22  A —  c2£>nF44 AAT^q  —  2ci  c2 Fi  1  F44A 

N xq  —  c2DqqF55XNx0  +  4a2ci£>i2F06 AAFco  _  4cic2D55FqqXNxo  +  8a2c2Fi2F66 AA^^O  _  2cic2F55F66A 

NxO  —  c*2c2£>22FJii  XNxq  —  2c2c2F44FFh  XNxq  —  2cx2  c\f22HixXNxq  —  c2c2F44FFn  XNxq  +  2a2c2Fi2FFi2  A 

A^ccO  +  ^cx2c\DqqHi2  XNxq  +  4a2cf  Fi2FFi2  AAT'jco  +  4a2cf  FqqHi2  XNxq  +  a2  c4FF22  A  AT^o  —  a2c2£>nFJ22  A 

^ccO  —  2a2cf  FhFF22AA/'£Co  —  a2c4FFn  FF22  AAT^q  +  2a2c2Fi2FF66  AAT^q  —  2clc2D55H66XNx0  +  4a2c^Fi2FF66 

AAFcO  -  clc2F55H66XNx0  +  2cx2 cl Hi2 Hqq XNxq  -  A44  (a55Fh  +  2c2FnD55  +  c2£>iiF55  +  a2 (2ci  (F±1  (D12  +  2F66 
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—  Cl  (F12  +  2 F66))  +  ci(Di2  +  2Dqq)H11)  +  Dh(Dqq  +  ci(-2(Fi2  +  Fqq)  +  ciHqq)))  +  (-Du  +  Cl  (2Fn  +  Ci  ifn  ))  A-ZV-eq) 
+  ^55  (  —  C2D11  (2D44  +  C2F44)  +  a2  (D22  —  Dh(D22  +  ci(2F22  +  ci  H22))  +  2D12 (Dee  +  ci(Fi2  —  c\Hqq)) 

+  ci  (ci(Fi2  +  2 F66)2  +  4D66(Fi2  -  ci  if 66))  )  -  (D66  +  ci  (2F66  +  c\ Hqq))  AA7'£c0)  -  ck2(Dh  +  ci(2Fn  +  ciifu)) 

(De6  +  ci(2F66  +  ciHqq))XNvq 

a4  =  4A55a2C2-Di2-D44  —  2a4C2D22D44  +  2a4C2DiiD22D44  +  8a2C2Di2D44D55  +  A55a4  D22Dqq  +  8A55a2C2D44D66 

—  4a4C2Di2D44DQ6  +  2a4  c2D22DqqDqq  +  16a2C2D44D55De6  +  4a4ciC2D22^44Fn  —  2A^,r)Oi4c\D22Fi2  —  4a4c\c2D\2 
D44F12  —  4a4ciC2D22D55Fi2  —  8a4c\c2D44DQQF\2  —  a® c\d22F^2  ~  2a4  c\c2D44Ff2  +  2A55 a4ci D12 F22  +  4a4ci 
C2D12D55F22  +  4A55oi4ciD66F22  +  8a4cic2D55D66F22  —  2A55a4c2Fi2F22  +  2a6c2Di2Fi2F22  —  4a4c2c2D55Fi2 

F22  +  2a6c?F22F22  -  ot&c\DlxFl2  -  2o?c\f11F%2  +  2A55a2  cId12F44  -  ol4 c\d\^F4 4  +  a4^!!^^  +  4a2 
c2  D12  D55  F 44  +  4A55a2  c2DqqF44  —  2a4  c2D12DqqF44  +  8a2C2D55D66F44  +  2a4ciC2l>22^1lf,44  —  2a4ciC2Di2 
F12F44  —  4a4  ci  c^DqqF\2F44  —  a4c\c^F22F44  +  4a2c2Di2D44F55  +  a4  c^D^DqqFqq  +  8a2c2D44D66D55  —  2a4 
cic2D22F12F55  +  2a4c1C2D12F22F55  +  4a4cic2D66f22i?55  —  2a4  c2c2F12F22Fqq  +  2a2C2Di2F44F55  +  4a2 
c2D66F44F55  —  2A55a4  ciD22Fqq  —  4a4  cic2D22D55F66  —  4a6 c2  D22F12Fqq  —  8a4 c2c2D44F12  Fqq  —  4A55a4c2F22 
Fqq  +  4a6c2Di2F22D66  —  8a4  c2c2D55F22Fqq  —  4a4  c2C2F12F44Fqq  —  2a4  c\c^D22FqqFqq  —  4a4c\c^F22FQQ 
F 66  -  4a6c2D22D26  -  8a4 c\c2D 44F^  -  8aQ c\f22F2q  -  4a4c2c2F44D26  +  2  a4  c\c2D22D  44H1X  -  a6c4 
F22  H 1  4  +  a4c\c^D22F44H\\  +  2a6  c2D22DqqHi2  —  2a6c\D22Fi2Hi2  +  2a6cf  Di2F22Hi2  +  4aQ  c\DqqF22H\2 
+  2a 6  c4  Fi2F22  if  12  +  4a6c4F22D66iii2  —  cx6  c!d22H22  +  2  A55  a4c2Di2 -H22  —  a6  c2D22H22  +  a6c2Dn  D22H22 
+  4a4c2C2Di2D55-Ff22  +  4A55a4  c2DqqH22  —  2a6  c2D12DqqH22  +  8a4c2C2D55D66if22  +  2a6cf  D22FnH22  —  2  a6cf 
F)12F12H22  —  4a6  c\DqqF12H22  —  a6  c4F22H22  +  2a4  c2C2Di2Fq5H22  +  4a4  c\c^DqqFqqH22  —  4a6  c1f12FqqH22 

—  4a& c\FqqH22  +  a6 c4 D 22 H \  \  H 22  +  A^^a4  c2D22Hqq  +  4a4c2C2Di2D44-Ffe6  +  2a4c2c2D22DQQHQQ  +  4a6c2D22 
d66h66  +  8a4c2c2  D44DqqHqq  —  4a6  c^D22F12Hqq  +  4a6cf  D12F22Hqq  +  8a6cf  DqqF22Hqq  +  2a4c2C2Di2F44ff66 
+  4a4c2C2D66F44-ff66  +  a4c2c2D22F^H6Q  —  2a6  c1d22H12Hqq  +  2a6  c!d12H22Hqq  +  4a6c4D66ff22ff66  +  ^55«2 
D22XNx0  +  2A55C2D44AiVxo  +  2a2  c2D22D55XNxq  +  4C2D44D55  XNx0  +  a4  D22Dqq  XNxq  +  2a2  c2D44DqqXNx0 

+  2A55a2ciF22  A-ZV^o  +  4a2ciC2D55F22  A-ZV^q  +  2a4c\DQQF22  A/V^o  +  ^55c2  i744^-^x0  +  2  c2  D  5  5  F44  A  iV^  0  +  a2 
02-^66-^44  A-ZV-^o  +  ol2c^D22F5$XNxq  +  2C2D44F55  XNxq  +  2a2  ci  c2  F22  D55  A  +  C2F44F55  XNxq  +  2a4ci 

D22  FqqXNxO  +  4a2ci  c2  D  44FqqXN  xq  +  4a4c2F22D66AA7'a,o  +  2a2  c\c^F44FqqXNx0  +  A55a2c2H22  XNx0 
+  2a2c2c2D55H22XNx0  +  a4  c2D6qH22XNx0  +  a2  c2C2F55H22XNx0  +  2a4cf  FqqH22  XNx0  +  a4c2D22 
HqqXNxO  +  2a2c2C2  D44HqqXNxq  +  2a4  c\f22HqqXNxq  +  a2c\c^F44HQQ  XNx0  +  a4c4ii22if66^-^a;0  +  a2 
(^4-55 (-£>66  +  ci(2F66  +  CiHqq))  +  2c2(Dn  D44  +  DqqDqq  +  ci(2D44Fii  +  2DqqFqq  +  C1D44F11  +  c\DqqHqq)) 

+  C2(DhF44  +  DqqFqq  +  ci  (2Fn F44  +  2 F55F66  +  ci  F44 H !  4  +  ciF55if66))  -  a2  (d22  -  £>n(D22  +  c1(2F22 
+  ci  H22))  +  2D12(Dqq  +  ci(2Fia  +  2Fqq  +  ci  (if  12  +  if  66)))  +  ci(4Fi2(D66  +  ciFi2)  -  D22(2Fn  +  ciifu) 

+  ci(8Fi2F66  +  2DqqH\2  —  2Fn  (2F22  +  c\H22)  +  ci(  —  2F22Hn  +  ifi2(4(Fi2  +  F66)  +  ciiii2)  —  ciHhH22 
+  4Fi2 if66  +  2ciiii2if66)))))^ify0  +  ^44  (2A55a2(Di2  +  2Dqq)  +  a4  (  —  D22  +  Di  1  D22  +  ci  ^22^11  —  4D66 
i^l2  -  ci(Fi2  +  2Fqq)2  +  ci D22Hh  +  4ci DqqHqq^  -  2D12(Dqq  +  ci(F12  -  ci if66)))  +  A55 A-ZV-eq  +  c2(2D55  +  c2 
F55) XNxo  +  a2  ^4c2 D55 (D12  +  2Dqq)  +  2c2(Di2  +  2Dqq)Fqq  +  (D66  +  ci(2F66  +  ciHqq))XNxq  +  (Du  +  ci 
(2Fn  +  ci  ifii ))  A-ZVyo)  ) 


( 


(^55  +  2c2D55  +  C2F55  +  a2  (d66  +  2c  1  Fqq  +  c2ifee))  (a2  (^44if22  +  2C2D22D44  +  C2D22D44  +  a2c2 
+  D22H22fj  )  +  (a44  +  c2  (2D44  +  C2F44)  +  a2  ^D22  +  2ci  F 22  +  c2if22))  ^ffyo) 


(B.l) 
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Appendix  C:  Gantt  Chart 
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Figure  11:  Gantt  Chart 
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